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LETTER  FROM  THE  EDITOR-IN  CHIEP.  ACES  JOURNAL 


We  live  in  times  when  our  personal  and  professional  circumstances  are  increasingly  subject  to  rapid 
change.  The  same  is  true  of  the  ACES  Journal.  The  March  1994  issue  is  the  first  one  which  I  have  had 
the  privilege  of  putting  together  as  the  new  and  untried  Editor-In-Chief.  I  now  realize  the  tremendous 
contributions  that  David  Stein  has  made  over  the  years  as  Editor-In-Chief  for  ACE^.  This  contribution 
can  be  measured  against  the  fact  that  it  has  taken  three  people,  Perry  Wheless,  Jr.,  Adalbert  Konrad  and 
myself  to  replace  him.  Perry  is  the  new  Editor-in-Chlef  for  ACES  with  overall  responsibility  for  the  ACES 
publications,  while  Adalbert  is  Associate  EdItor-in-Chief  of  the  ACE^  Journal.  Adalbert's  expertise  will 
be  particularly  useful  in  the  field  of  low  frequency  power  electromagnetics.  On  behalf  of  our  Joumzd  1 
acknowledge  and  thank  Dave  Stein  for  his  dedication  and  hard  work  over  the  years  to  ensure  that  the 
Journal  will  occupy  a  position  of  pre-eminence  in  the  applied  electromagnetics  commuirity. 

This  year  also  marks  the  Tenth  Review  of  Progress  in  Spiled  Computational  Electromagnetics.  It  is 
appropriate  that  ACES  and  all  its  committees  should  take  stock  of  their  achievements  and  plan  a  new 
course  for  the  future  at  this  time.  The  Journal  is  always  reviewing  its  role  within  ACE^  and  investigating 
ways  and  means  of  better  serving  its  readership.  With  the  able  assistance  of  the  panel  of  Journal  editors. 
Perry.  Adalbert  and  I  are  reviewing  the  type  of  article  which  would  have  the  greatest  value  for  our  readers. 
Inevitably,  there  will  be  differences  of  opinion  because  of  the  diversity  of  interests  of  the  ACE^S  members 
and  readers  of  the  Journal.  We  will  strive  for  the  best  possible  compromise  to  ensure  that  the  articles  will 
be  of  Interest  to  the  broadest  cross-section  of  our  members  and  readers,  and  that  these  articles  will  support 
the  objectives  of  ACES.  As  Editors-in-Chlef  we  do  however,  reserve  the  right  to  publish  articles  which  may 
be  contentious  and,  hopefully,  thought  provoking  as  well  as  those  which  may  open  up  a  new  area  of 
investigation. 

As  with  all  technical  publications,  the  ACEiS  Journal  and  Newsletter  are  faced  with  ever  increasing 
publication  costs.  This  issue  serves  as  a  trial  run  for  one  method  of  reducing  these  costs.  The  reader  will 
notice  that  several  articles  have  been  published  in  a  compact  format  using  double  columns  and  a  reduced 
font  size.  The  effect  of  this  has  been  to  reduce  the  number  of  printed  pages,  and  thus  costs,  without 
affecting  the  actual  length  of  the  papers  themselves.  Our  sincere  thanks  to  those  authors  who  have  been 
able  to  assist  with  this  trial  run.  Our  readers  are  invited  to  give  constructive  criticism  of  these  changes. 
Please  note  that  no  author's  paper  will  be  turned  down  because  he/ she  is  unable  to  accomodate  the 
compact  format. 

As  Editor-in-Chlef  of  the  Journal  I  gratefully  acknowledge  the  support  and  assistance  of  Perry  Wheless, 
Adalbert  Konrad,  Dick  Adler  and  Dave  Stein  in  putting  this  issue  together.  Not  surprisingly,  Dave  has  been 
particularly  generous  in  passing  on  the  wisdom  and  experience  he  has  built  up  over  the  years. 

Finally,  the  Journal  reflects  the  collective  wisdom  and  hard  work  of  each  and  every  author,  along  with  the 
constructive  criticism  of  various  editors  and  unnamed  reviewers  who  have  given  generousty  of  their  time 
and  experience.  Without  them  there  would  be  no  Journal. 

Duncan  C.  Baker 
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ON  THE  USE  OF  BIVARIATE  SPLINE  INTERPOLATION  OF  SLOT  DATA 
IN  THE  DESIGN  OF  SLOTTED  WAVEGUIDE  ARRAYS 

Derek  A.  McNamara  and  Johan  Joubert 

Depanment  of  EUarical  <&  Electronic  Engineering.  University  of  Pretoria, 

Pretoria,  South  ^nca  0002 


Abstract  :  Essential  to  the  design  process  of  slotted 
waveguide  arrays  is  an  accurate  knowledge  of  the  self- 
properties  of  the  individual  slots.  These  properties  can  be 
computed  with  sufficient  accuracy  using  moment  method 
codes,  which  are  too  computationally  intensive  to  be 
incorporated  into  design  codes  for  large  arrays.  This 
paper  describes  an  interpolation  scheme  which  provides 
both  sufficiently  accurate,  and  rapid,  computation  of  slot 
data  from  a  pre-determined  database  to  be  viable  for  use 
in  a  slotted  waveguide  array  design  code. 

1.  INTRODUCTION 

It  is  correct  to  state  that  the  structured  design  of 
high-performance  slotted  waveguide  arrays  began  with 
the  publication  by  Elliott  [1]  in  1978.  All  the  design 
procedures  available  in  the  literature  are,  to  the  best  of 
the  authors’  knowledge,  based  on  this  original  approach 
(2,3],  with  further  discussions  [3]  aimed  principally  at 
examining  the  correlation  of  the  design  method  with 
experiment  or  the  details  of  its  implementation. 

Essential  to  the  design  process  is  an  ability  to 
determine  the  self-properties  of  the  individual  slots. 
These  properties  can  be  computed  with  sufficient 
accuracy  using  moment  method  techniques.  In  a  recent 
paper  [3],  which  very  carefully  compares  theory  and 
experiment,  the  inclusion  of  the  moment  method 
computation  itself  in  the  design  process  has  been 
suggested.  However,  because  of  the  extensive 
computational  effort  involved,  this  has  only  been  used  in 
connection  with  the  design  of  linear  arrays  of  slots  with 
no  more  than  21  elements.  An  even  braver  approach  has 
been  presented  by  Gulick  and  Elliott  [4],  where  the 
moment  method  is  used  to  obtain  a  complete  solution  for 
the  integral  equation  formulation  of  the  array  antenna 
considered  as  a  whole.  The  advantage  of  this  scheme  is 
that  it  does  not  require  individual  characterisation  of  the 
radiating  elements,  with  mutual  coupling,  both  internal 
and  external,  being  handled  automatically.  However,  this 
approach  has  thus  far  only  been  demonstrated  on  an  array 
of  8  slots  and  is  at  present  not  a  practical  procedure  for 
the  design  (as  opposed  to  analysis  only)  of  large  planar 
arrays. 


Although  the  moment  method  formulations 
(S,6,7]  are  quite  efficient,  they  are  still  too 
computationally  intensive  to  be  incorporated  into  design 
codes  for  large  arrays.  Although  it  is  possible  to  use 
approximate  formulas  for  slot  properties,  it  is  possible  to 
avoid  any  compromise  and  use  numerically  generated  slot 
data  along  with  a  spline  interpolation  scheme  for  use  in 
the  design  phase.  This  has  already  been  an  approach 
recommended  in  microstrip  circuit  design  [8]  and  the 
modeling  of  antenna  near  a  half-space  [11]. 

When  designing  slotted  waveguide  arrays,  the 
waveguide  dimensions  (e.g.  full-height,  half-height, 
quarter-height)  and  wall  thickness  would  have  been 
determined  before  the  start  of  the  design.  Practical 
fabrication  considerations  usually  require  that  the  slot 
width  be  the  same  pre-determined  value  for  each  of  the 
slots.  Also,  although  the  array  performance  may  be 
analysed  at  a  number  of  frequencies,  any  design  run  is 
done  at  some  single  frequency.  Thus,  if  the  results  of  the 
moment  method  code,  namely  complex  s,,  and  Sji  (or 
complex  self-immittances  of  the  equivalent  network 
model)  are  obtained  over  the  range  of  slot  offsets  x,  and 
lengths  2t  (these  are  shown  in  Fig.l),  such  sets  of  data 
can  be  considered  as  bivariate  functions  of  the  variables 
Xg  and  21.  Ready-to-use  computer  routines  are  available 
[9]  for  performing  a  bivariate  spline  interpolation,  and 
these  have  been  used  by  the  present  authors. 

2.  BIVARIATE  CUBIC  SPLINE  INTERPOLATION 
OF  SLOT  SCATTERING  PARAMETERS 

The  procedure  is  to  establish  (via  the  moment  method 
code)  the  slot  data  at  some  single  frequency  at  selected 
values  of  the  two  independent  variables  x„  and  2t .  The 
range  of  this  data  will  usually  be  known  from  previous 
experience  or  examination  of  the  literature.  This  data 
may  be  the  complex  Sn  and  S21,  from  which  equivalent 
network  parameters  are  easily  obtained,  or  such 
equivalent  network  parameters  themselves.  The  routine 
SPUE2  [9,P.  100]  is  then  run  with  this  data  as  input  in 
order  to  compute  certain  derivative  data  needed  for  the 
final  spline  interpolation.  Thereafter  the  original  slot 
data,  plus  that  generated  by  SPIJE2,  is  used  with  routine 
SPIJN2  [9,p.lOI]  to  compute  interpolated  values  of  the 
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slot  data  for  any  given  x,,  and  21  combination  as 
requested  by  a  design  code.  This  process  is  as  rapid  as  if 
direct  closed-form  formula  such  as  the  variational 
expression  in  [10]  were  utilised  (which  is  not  sufficiently 
accurate  of  course),  even  when  many  data  nodes  are 
used.  The  latter  remark  applies  to  the  other  quantities  as 
well.  Note  that  the  above-mentioned  routines  also  require 
the  availability  of  two  others  {SPLINE  and  SPLINT),  also 
available  in  [9,pp.86-89]. 

In  order  to  illustrate  this  point,  Table  1  shows 
the  scattering  parameters  (magnitude  and  phase)  for 
various  slot  lengths  and  offsets  (at  9.0  GHz  in  half-height 
X-band  waveguide)  computed  with  the  moment  method 
code  and  with  the  spline  interpolation  procedure  from  a 
database.  The  width  of  the  slots  were  taken  as  1.6  mm 
and  the  wall  thickness  of  the  waveguide  as  0.9  mm. 
Database  values  were  computed  with  the  moment  method 
code  at  23  node  points  in  slot  length  (13  mm  to  24  mm 
in  steps  of  .5  mm)  and  12  node  points  in  slot  offset  (.8 
mm  and  1  mm  to  6  mm  in  steps  of  .5  mm).  The  database 
was  computed  to  an  accuracy  of  6  significant  digits,  so 
as  not  to  lose  accuracy  during  the  interpolation  process 
through  round-off  error.  All  the  lengths  and  offsets 
shown  in  Table  1  are  in  between  these  node  points  and 
examination  of  the  data  reveals  that  the  difference 
between  the  exact  moment  method  computed  data  differ 
only  very  slightly  from  the  scattering-parameter  data 
computed  from  the  database.  Fig. 2  has  been  included  to 
give  a  global  view  of  the  variation  of  the  magnitude  of  s,, 
with  the  slot  parameters  x„  and  2i.  We  have  shown  a 
portion  of  the  x,-3mm  dataset  in  Fig. 3,  obtained  using 
a  decreasing  number  of  nodal  points  to  give  some 
indication  of  the  number  needed.  The  accuracy 
achievable  can  easily  be  increased  simply  by  using  more 
nodes  in  the  database  until  it  is  well  within  that  required 
for  design  purposes  (and  certainly  up  to  that  accuracy 
with  which  slot  properties  can  be  measured). 

3.  CONCLUDING  REMARKS 

The  interpolation  schemes  whose  use  has  been 
discussed  in  this  paper  provide  both  sufficiently  accurate, 
and  rapid,  computation  of  slot  data  for  their  use  in  a 
slotted  waveguide  array  design  code  to  be  viable.  Given 
a  set  of  measured  or  computed  (using  accurate  moment 
method  techniques)  data,  node  points  can  be  selected  in 
such  a  way  as  to  obtain  sufficient  accuracy  with  as  few 
nodes  as  possible  (and  thus  interpolate  as  rapidly  as 
possible). 
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Data  computed  with  the  moment 
method  code 

Data  con 
spline  ini 

iputed  froB 
terpolation 

a  database 

using 

2t 

B 

|Snl 

Zs,, 

|s„l 

ZSr. 

|S„I 

Z^it 

|S:.| 

ZS:. 

13.2 

2.1 

.027 

76.05 

.994 

-1.62 

.027 

76.00 

.994 

-1.63 

14.4 

3.4 

.104 

63.76 

.959 

-5.93 

.104 

63.76 

.959 

-5.93 

15.6 

EB 

.188 

43.46 

.874 

-9.21 

.188 

43.46 

.874 

-9.21 

16.8 

4.1 

.284 

15.47 

.732 

-7.37 

.284 

15.67 

.733 

-7.43 

17.2 

4.8 

.337 

15.41 

.684 

-9.85 

.337 

15.51 

.685 

-9.87 

18.7 

5.2 

.380 

-8.46 

.622 

0.74 

.380 

-8.44 

.622 

0.73 

Table  1:  Comparison  between  scattering  parameter  data  generated  by  a  moment  method  code,  and 

data  computed  firom  a  data  base  using  a  spline  interpolation  routine.  (Lengths  and  offsets 
given  in  mm,  and  phase  in  degrees). 


Figure  1:  Geometry  of  a  longitudinal  slot  in  rectangular  waveguide,  showing  the  notation  used  in 

the  paper. 
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Magnitude  of  sll 


Figure  2;  Magnitude  of  s,,,  for  varying  x„  and  2f,  computed  using  the  moment  method. 


Figure  3:  Magnitude  of  s,,,  for  x,=3nun  and  varying  2t,  computed  using  the  spline  interpolation 

procedure  with  different  numbers  of  nodal  points.  Some  moment  mediod  results  are 
shown  (*)  for  reference. 
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ABSl  RACT:  A  study  of  the  determination  of  non- integer  eigenvalues  for  ordinary  differential 
equations  with  transcendental  solutions  is  presented.  An  algorithm  based  on  expansion  in  terms 
of  ( 'hebyshev  polynomials  and  collocation  is  presented.  The  method  is  applied  to  the  problem 
of  computing  the  electric  field  external  to  a  biconical  radiating  structure.  Eigenfunction  solutions 
for  Legendre  s  differential  equation  satisfying  the  boundary  conditions  of  the  problem  considered 
are  presented. 


1.  INTRODUCTION 

Canonical  analysis  of  boundary  value  problems  commonly  produce  series  solutions  of 
transcendental  functions,  fl-3]  where  the  index  of  summation  is  a  set  of  eigenvalues  (v) 
determined  by  the  boundary  conditions  of  the  problem.  When  considering  perfectly  conducting 
boundaries,  the  eigenvalues  are  those  values  such  that  the  transcendental  function  of  order  v,  or 
its  derivative,  is  equal  to  zero  at  the  boundaries.  Unfortunately,  the  determination  of  the 
eigenvalues  that  provide  such  results  is  usually  limited  to  special  cases,  often  when  v  is  an 
integer.  Difficulty  in  determining  the  correct  values  for  v,  when  v  is  not  an  integer,  has  limited 
the  use  of  canonical  analysis  for  these  problems. 

This  paper  presents  a  generalized  numerical  approach  for  determining  the  eigenvalues  of 
transcendental  functions  subject  to  the  Dirichlet  and  Neumann  boundary  conditions.  The  method 
proposed  is  based  upon  expanding  the  unknown  function  in  a  series  of  Chebyshev  polynomials 
(4]  and  using  the  method  of  collocation  [3]  to  obtain  a  well  conditioned  system  of  linear 
homogeneous  equations.  The  eigenvalues  (v)  are  then  found  by  using  a  bracketing  and  bisection 
technique  [5]. 
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This  paper  is  organized  as  follows.  Section  2  describes  the  mathematical  formulation  of 
the  problem.  Section  3  presents  the  method  for  numerically  computing  the  eigenvalues.  Section 
4  presents  sample  results  for  Legendre’s  differential  equation,  which  arises  in  the  canonical 
analysis  of  spherical  based  problems. 

2.  MATHEMATICAL  FORMULATION 

Assume  that  the  physical  problem  under  consideration  is  described  by  the  function  y„(x), 
which  is  defined  on  the  closed  region  [a,b],  and  satisfies  a  second  order  ordinaiy  differential 
equation  (ODE)  of  the  following  form: 

dy  (x) 

A(x)  +  Bix)  +  Cix,v)  y^(x)  =  0  .  (2.1) 

dx^  dx 

The  functions  A(x),  B(x),  and  C(x,v)  are  taken  to  be  continuous  on  the  open  region  (a,b).  The 
sought  solutions  are  subject  to  boundary  conditions  which  can  be  separated  into  different  types 
according  to  the  physical  problem  under  consideration.  Typical  boundary  conditions  are: 


yv(a)  =  yv(b)  =  0  (2.2) 

y;(a)  =  y;(b)  =  0  (2.3) 

yv(a)  =  y;(b)  =  0  (2.4) 

y;(a)  =  y,(b)  =  0  (2.5) 

where  yv'(^)  denotes  the  derivative  of  yv(x)  with  respect  to  x,  evaluated  at  x  =  ^. 

Let  y^,(x)  be  expressed  by  the  following  Chebyshev  polynomial  expansion: 

00 

yvW  =  E  ««  T„[(!(a:)]  a^x^b  (2.6) 

n=0 


where  T„(z)  is  the  nth  Chebyshev  polynomial  of  the  first  kind,  and  f  (x)  is  a  linear  mapping  which 
maps  the  interval  [a,b]  to  [-1,1] 


II 

s 

cos(/i  cos  '(z)) 

(2.7) 

f(A:)  = 

2  2a  ^ 

(b-a)  ^  (b-a) 

a^x^b 

(2.8) 

Since  the  set  of  Chebyshev  polynomials  are  continuous  on  [-1,1],  the  first  and  second  derivatives 
of  y^,(x)  may  be  expressed  as  follows: 
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(2.9) 


dx 


E  <  <  [«w] 


d((x) 

dx 


d\ix) 

dx^ 


=  E  “»  T"  [«(j:)j 


(2.10) 


T|^[J(j:)]  and  T^^[{(a:)]  are  the  first  and  second  derivatives  of  T„[{(jc)]  with  respect  to  f  (x),  and  are 
given  by: 

Tte  =  — [-nz  T„(z)  +  n  T,.,(z)]  (2.11) 

(1-z  ) 

[  T„(z)  {(nz)^  -  nz^-n) 

(1-zT 

+  T„.,(z)(-2n2z  +  3nz)  +  T„.2(z)(n2-/i)  ]  .  (2.12) 

Substituting  Eqs.  (2.9)  and  (2.10)  in  Eq.  (2.1)  and  approximating  the  series  expansion  for  y^(x) 
by  the  first  N  terms,  yields  the  following  linear  homogeneous  equation  for  the  expansion 
coefficients  (a^): 


E 


<(x)  a;  =  0 


(2.13) 


where, 


M:(x)=A(x)T''[llix)] 


dm 

dx 


2 

+ 


5(x)T;[(!(x)] 


dm 

dx 


+  C(x.v)T„[f(x)] 


(2.14) 


Enforcing  Eq.  (2.14)  at  N  points,  {xj,  i=l,  N},  on  the  interval  [a,b]  leads  to  a  system  of  N  linear 
homogeneous  equations,  which  may  be  written  as  the  following  matrix  equation: 
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•  •  *  ^W-l('^l) 

V 

«o 

0 

•  •  • 

V 

“l 

— 

0 

•  •  • 

V 

“n-1 

0 

(2.15) 


Note  that  the  matrix  M  is  solely  a  function  of  v,  the  desired  eigenvalues  of  the  given 
ODE.  Since  Eq.  2.13  is  a  homogeneous  equation,  it  will  have  non-trivial  solutions  if  and  only 
if  the  determinant  of  the  matrix  M  is  zero.  However,  taking  the  determinant  of  the  matrix  M  in 
its  present  form  yields  zero  for  any  value  of  v.  This  is  because  solutions  exist  for  any  given 
value  of  V  due  to  the  fact  that  the  boundary  conditions  have  not  yet  been  imposed.  Hence,  it  is 
necessary  to  first  impose  the  boundary  conditions  in  order  to  obtain  the  desired  values  for  v. 

Boundary  conditions  are  imposed  by  replacing  the  first  and  last  rows  of  the  matrix  3/ with 
the  series  representation  for  the  boundary  conditions.  Thus  if,  y„(a)  =  y^Cb)  =  0  the  first  row  of 
the  matrix  M  is  replaced  by 

E  <  =  0  (2.16) 

«=0 

and  the  last  row  is  replaced  by 

E  T,[((i.)]  =  0  ,  (2.17) 

n=0 

The  new  matrix  obtained  will  be  denoted  by  M.  Because  the  boundary  conditions  require  the 
function,  or  the  derivative  of  the  function,  to  be  zero  at  the  boundaries,  the  matrix  equation 
remains  homogeneous;  thus,  non-trivial  solutions  still  exist  if  and  only  if  the  determinant  of  the 

matrix  M  is  zero.  Hence,  the  permissible  values  of  v,  subject  to  the  given  boundary  conditions, 
are  obtained  by  requiring  det(Af)  =  0. 
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3.  COMPUTER  IMPLEMENTATION 


Numerical  estimation  of  the  eigenvalues  (p  =  l,...,P),  is  based  on  the  fact  that  the 
det(iW)  is  an  oscillatory  function  of  v,  with  the  det(M'^'’)  =  0  (p  =  Hence,  over  any 

interval  containing  a  single  eigenvalue,  det(M)  will  change  sign.  This  behavior  allows  the 
eigenvalues  (vp  to  be  determined  using  a  bracketing  and  bisection  technique  [5].  Scanning 
det(Af)  for  changes  in  sign  over  a  given  interval  on  the  v  axis,  provides  the  bracketing  intervals 
for  the  eigenvalues.  Care  must  be  taken  in  selecting  a  maximum  scan  distance  which  is  less  than 
the  minimum  distance  between  any  two  adjacent  roots.  Scan  distances  which  are  too  large  may 
cause  roots  to  be  missed.  Once  the  roots  are  bracketed  the  bisection  technique  may  be 
implemented  to  compute  the  particular  eigenvalue  to  the  desired  degree  of  accuracy.  Since  the 
bisection  method  requires  only  the  computation  of  the  sign  of  the  determinant,  the  common 
problem  of  numerical  over-flow,  associated  with  determinant  calculations,  is  avoided. 

4.  ANALYSIS  OF  A  BICONICAL  RADIATING  STRUCTURE 

The  technique  developed  in  section  2.  is  now  applied  to  Legendre’s  differential  equation 
which  describes  the  electromagnetic  field  external  to  a  biconical  radiating  structure  (shown  in 
figure  4. 1 ).  It  can  be  shown  that,  under  radiation  conditions,  the  electric  fields  external  to  the 
structure  have  series  solutions  of  the  following  form: 

=  E  -ffvM/2(^)  ^v(COS0)  (4.1) 

V 

(2) 

where  //y^,/2(kr)  is  a  modified  Hankel  function  of  the  second  kind,  and  L„(cos0)  is  an  odd 

Legendre  polynomial  [6].  Boundary  conditions  require  that  (cos0,)  =  Z,„  (cos02)  =  0.  Thus, 
it  is  necessary  to  determine  the  eigenvalues  whose  eigenfunctions  satisfy  the  given  boundary 
conditions.  For  Legendre’s  differential  equation,  A(x)  =  l-x\  B(x)  =  -2x  and  C(x,v)  =  v(v+l). 

Two  different  biconical  radiating  structures  were  chosen  for  analysis:  Structure  1,0,= 
30°  and  02  =  150°;  Structure  2,  0,  =  10°  and  02  =  170°.  Tables  I  and  II  provide  the  first  four 
eigenvalues  for  each  structure  emd  compares  the  calculated  eigenvalues  with  the  eigenvalues 
estimated  by  Grimes  using  an  asymptotic  expansion  technique  [7].  Graphs  of  the  corresponding 
eigenfunctions  are  shown  in  figures  4.2  and  4.3.  For  both  radiating  structures  under  study, 
coverage  of  the  calculated  eigenvalues  occurred  for  matrix  sizes  on  the  order  of  N  =  50  to  60. 
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Figure  4.1  Geometry  of  the  biconical  radiating  structure 


TABLE  1 

Structure  1  =  30%  6^  =  ISO"* 


Eigenvalue 

Calculated  Value 

Grimes’  Result 

nu-1 

2.439211 

2.439211 

nu-2 

5.466996 

5.466996 

nu-3 

8.477510 

8.477309 

nu-4 

11.482784 

- 

TABLE  n 

Structure  2  (?,  =  10°,  =  170° 


Eigenvalue 

Calculated  Value 

Grimes’  Result 

nu-1 

1.621407 

1.620624 

nu-2 

3.916836 

3.915488 

nu-3 

6.188799 

6.187171 

nu-4 

8.451585 

8.450112 

LEGENDRE  FUNCTIONS  «  LEGENDRE  FUNCTIONS 


Figure  4.3  First  four  Legendre  functions  (eigenfunctions)  for  30®  biconical  structure 

(0,  =  30®,  02  =  150®) 
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5.  CONCLUSIONS 


A  generalized  technique  for  determining  non-integer  eigenvalues  for  ordinary  differential 
equations  with  transcendental  solutions  has  been  presented.  To  verify  the  technique,  the  method 
was  applied  to  a  biconical  radiating  structure  and  the  computed  eigenvalues  were  compared  with 
results  obtained  by  Grimes  using  an  asymptotic  method.  Excellent  agreement  was  established 
between  the  two  methods.  The  method  proposed  here  is  completely  general  and  has  the 
advantages  that  only  the  sign  of  the  determinant  needs  to  be  computed  and  that  the  matrix  size 
needed  for  convergence  is  small. 
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Abstract 

This  paper  describes  the  electrical  characterization  of  the  Very  Low  Frequency  (VLF) 
dual  trailing  wire  antenna  systems  associated  with  Navy  aircraft  used  for  communications  to 
underwater  Navy  assets.  The  static  wire  geometry  with  the  aircraft  in  orbit  was  obtained 
using  a  modified  version  of  the  program  WIRES,  which  models  the  mechanical  parameters  of 
the  dual  trailing  wire  antenna.  The  Numerical  Electromagnetics  Code  (NEC)  was  used  to 
model  all  electrical  characteristics  of  the  VLF  dual  trailing  wire  systems  associated  with  the 
typical  TACAMO  aircraft.  Antenna  parameters  such  as  current  distribution,  input  impedance 
versus  frequency,  effective  radiated  power,  and  the  effect  of  antenna  wire  conductivity  were 
evaluated.  It  is  demonstrated  that  the  computer  codes  provide  results  which  are  consistent 
with  measured  values  of  static  antenna  electrical  parameters. 

I.  Introduction 

The  Numerical  Electromagnetics  Code  (NEC)  was  utilized  to  computer  model  the 
electrical  characteristics  of  the  VLF  dual  trailing  wire  transmitting  antenna  system,  which 
consists  of  both  a  long  and  a  short  trailing  wire  [1].  A  drogue  is  attached  to  the  far  end  of 
the  LTWA  which  causes  the  wire  to  be  as  vertical  as  possible.  A  drogue  is  also  attached  to 
the  end  of  the  STWA  in  order  to  provide  stability.  A  measure,  called  verticality,  is  defined 
as  the  ratio  of  the  projection  of  the  wire  on  the  Z-axis  to  the  total  length  of  the  wire. 
Verticality  is  given  in  terms  of  a  percentage,  and  the  higher  the  number,  the  better  for 
coupling  into  the  earth-ionosphere  waveguide  for  propagation.  The  wire  shape  geometries  of 
the  long  trailing  wire  antenna  (LTWA)  and  the  short  trailing  wire  antenna  (STWA)  were 
determined  using  the  programs  W1RE4NEC  and  STWANEC,  respectively  [2].  These  are 
steady-state  mechanical  codes  which  provide  piecewise  wire  segments  for  data  input  to  NEC. 
The  following  antenna  characteristics  have  bwn  obtained  for  this  system:  current 
distribution,  input  impedance  versus  frequency,  and  efficiency.  The  computed  antenna  input 
impedance  is  compared  to  measured  values  for  eight  actual  flights. 
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n.  Static  Characterization  of  the  Dual  Trailing  Wire  Antenna 

A  steady-state  mathematical  model,  which  computes  the  spatial  configuration  and 
tension  along  an  orbiting  LTWA,  was  originally  developed  at  the  Naval  Air  Development 
Center  (NADC),  Warminster,  PA  [3].  This  model  involved  the  numerical  solution  of  a 
system  of  three  coupled  ordinary  differential  equations  which  describe  the  position  of  the 
LTWA  with  respect  to  a  cylindrical  coordinate  system.  The  numerical  technique  used  to 
evaluate  this  system  of  differential  equations  was  later  improved  by  others  at  NADC  and  the 
improved  version  of  their  code  was  given  the  name  WIRE3  [4].  A  corrected  and  modified 
version  of  the  WIRE3  code  was  then  created  for  use  with  NEC  in  the  work  reported  here  and 
was  given  the  name  WIRE4NEC.  The  WIRE4NEC  code  writes  a  file  containing  a  set  of 
NEC  compatible  geometry  inputs. 

A  brief  summary  of  input  parameters  necessary  to  run  WIRE4NEC  follows.  All 
input  parameters  are  contained  in  a  data  file  called  WIREIN.DAT: 

Line  1  Number  of  runs  to  be  executed  (integer  number) 

Line  2  Run  title  (20  characters) 

Line  3  True  airspeed  (kts),  aircraft  orbit  radius  (ft),  aircraft  altitude  (ft),  number  of 
wire  increments  (integer  number),  step  size  (integer  number),  output  flag  (0 
for  trouble  shooting,  1  for  formatted  output). 

Line  4  Estimated  separation  between  the  aircraft  and  the  drogue  (ft),  range  of 
estimated  separation  (ft). 

Line  5  Drogue  weight  (lbs),  drogue  diameter  (ft),  coefficient  of  drag  of  the  drogue. 

Line  6  Wire  density  (Ibs/ft),  wire  diameter  (ft),  wire  length  (ft),  coefficient  of 
friction  of  the  wire,  coefficient  of  drag  of  the  wire. 

Lines  2-6  are  repeated  for  multiple  runs. 

The  WIRE4NEC  program  evaluates  the  wire  model  at  five  different  separation  cases 
between  the  aircraft  and  the  drogue.  The  program  will  stop  iterating  at  each  separation  case 
after  31  iterations  or  if  the  solution  converges  so  that  the  antenna  feed  endpoint,  known  as 
LTWA  distance,  is  within  100  feet  of  the  aircraft.  Generally,  it  requires  die  full  number  of 
iterations  for  the  program  to  converge  to  a  solution.  In  order  to  test  the  sensitivity  of  NEC 
to  the  endpoint  of  the  LTWA,  solutions  less  than  1(X)  feet  were  accepted. 

Several  models  were  considered  for  the  short  trailing  wire  antenna  (STWA).  The 
first  model  treated  the  STWA  as  remaining  in  the  plane  of  the  flight  path.  Two  cases  for 
this  assumption  were  analyzed:  STWA  tangential  to  the  radius  of  the  flight  path,  and  the 
STWA  along  the  circumference  traced  by  the  flight.  Both  models  provided  only  marginal 
results.  Next,  the  wire  was  allowed  to  drop  at  a  constant  angle  with  respect  to  horizontal. 
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resulting  in  an  improvement  over  the  previous  models.  Finally,  a  code  called  STWANEC 
was  developed  which  determines  the  steady-state  geometry  of  the  STWA  for  input  to  NEC. 
The  STWANEC  program  numerically  solves  a  system  of  coupled  ordinary  differential 
equations  which  represent  the  mechanical  forces  on  the  antenna.  A  fourth-order  Runge-Kutta 
technique  was  used.  A  complete  description  of  the  mechanical  model  for  the  STWA  is 
beyond  the  scope  of  this  work  [2]. 

m.  NEC  Modeling  and  Results 

Determining  the  geometry  of  the  Long  Trailing  Wire  Antenna  (LTWA)  and  the  Short 
Trailing  Wire  Antenna  (STWA)  of  the  dual  trailing  wire  system  is  done  using  the  two 
aforementioned  programs  WIRE4NEC  and  STWANEC.  Both  of  the  programs  utilize  various 
flight  data  in  order  to  predict  the  STWA  and  LTWA  positions.  Some  of  the  parameters  that 
are  used  in  both  programs  include:  true  air  speed,  aircraft  orbit  radius,  aircr^t  altitude, 
estimated  separation  between  aircraft  and  drogue,  and  drogue  weight.  Based  on  all  the 
parameters,  the  program  determines  the  verticality  of  the  wire  and  the  separation  of  the 
predicted  wire  endpoint  from  the  aircraft.  In  addition,  the  geometry  inputs  for  use  in  NEC 
are  produced.  By  merging  the  outputs  from  these  two  programs,  the  entire  geometry  of  the 
two  trailing  wires  is  obtained. 

Although  using  the  outputs  of  WIRE4NEC  and  STWANEC  yields  the  basic  geometry 
of  the  long  and  short  wires  in  a  NEC  format,  the  simple  concatenation  of  the  two  files  does 
not  produce  the  final  form  of  the  model.  Each  section  of  the  antenna  must  be  analyzed  in 
order  to  add  the  correct  amount  of  segment  length  tapering  near  the  feed  point.  It  is  also 
necessary  to  move  the  short  wire  so  that  it  is  connected  to  the  endpoint  of  the  long  wire/feed 
section  using  the  NEC  "move"  and  "translate"  input.  After  this  movement,  the  short  wire  is 
rotated  to  make  it  tangential  to  the  circumference  of  the  flight  path  using  the  NEC  "rotate" 
input. 


To  obtain  accurate  results  from  the  model,  the  wire  segment  lengths  near  the 
feedpoint  must  be  tapered.  An  algorithm  has  been  developed  to  calculate  the  adjacent 
segment  length  ratios  and  the  number  of  segments  for  both  the  long  and  short  wires.  The 
modeling  approach  is  to  choose  very  short  segments  near  the  critical  feedpoint  region  and 
then  to  smoothly  taper  outward  along  the  wires  from  this  feedpoint.  The  process  provides 
excellent  detail  at  the  feed  point  as  well  as  throughout  the  entire  antenna  structure.  Figure  1 
shows  the  top  view  and  Figure  2  shows  a  side  view  of  the  NEC  wire  segmentation  used  for 
the  dual  trailing  wire  antenna  shown  in  Figure  3.  The  tapering  algorithm  and  its 
implementations  will  now  be  described. 

NEC  has  an  option  for  tapering  either  the  length  or  radius  of  adjacent  segments  on  a 
straight  wire.  It  requires  as  an  optional  input  the  ratio  of  the  length  of  a  segment  to  the 
length  of  the  previous  segment  in  the  wire.  It  also  requires  the  radius  of  the  first  segment 
and  the  last  segment  if  stepped  radius  tapering  is  desired.  For  the  work  described  in  this 
paper  only  length  tapering  of  the  segments  was  required  which  will  be  discussed  in  the 
following. 
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In  the  method  used  in  NEC,  once  given  the  end  coordinates  of  the  wire,  the  number 
of  segments.  Ns,  and  the  adjacent  ratio  factor,  R^,  one  can  calculate  the  length  of  the  first 
segment,  S,,  as: 


where  L  is  the  total  length  of  the  straight  wire.  If  the  ratio  factor  equals  1 ,  then 

K  =  1. 


(1) 


(2) 


which  implies 


S,  =  L  /  Ns  (3) 

which  is  the  normal  case  of  the  wire  being  divided  up  into  a  total  of  Ns  equal  length 
segments. 

A  more  useful  and  desirable  approach  to  this  segment  tapering  problem  is  to  choose 
the  first  segment  length,  L,,  and  the  last  segment  length,  Lf, ,  and  with  these  parameters 

compute  the  number  of  segments,  Ns,  and  the  adjacent  segment  ratio  factor,  R^,  given  the 
total  length  of  wire,  L.  If  the  computed  R^  is  less  than  a  factor  of  2  or  some  other  modeling 
guideline  criterion,  then  this  R^  can  be  used. 


There  have  been  some  attempts  at  this  problem  [5]  using  an  iterative  solution  of  (1). 
This  method  is  time  consuming  and  unnecessary  because  there  exists  an  exact  analytical 
solution  for  this  problem  which  will  be  described.  The  last  segment,  Sff ,  is 


(4) 


This  can  be  rewritten  and  solved  for  N.  as 


N. 


log(5^  /5,) 

- — -  +  1 

log  R^ 


(5) 


Since  Ns  is  an  integer,  one  must  take  the  nearest  integer  of  this  expression  when 
programming.  Substituting  (4)  into  (1),  rearranging  and  solving  for  R^  gives 

(S«  -  L) 


(6) 
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A  FORTRAN  program  has  been  written  incorporating  the  above  equations  and  was  used  in 
the  modeling  discussed  in  this  paper. 

The  final  STWA  geometry  model  was  determined  using  the  steady-state  mechanical 
code  called  STWANEC  discussed  above.  This  model  resulted  in  a  curv^  drop  for  the 
STWA.  Furthermore,  the  wire  dropped  such  that  it  was  tangential  to  the  flight  path 
circumference.  The  final  form  for  the  model  can  be  seen  in  Figure  3  which  plots  the  LTWA 
and  STWA  geometries  from  different  aspect  angles.  Using  this  final  form,  excellent 
agreement  was  obtained  with  respect  to  the  measured  data.  A  comparison  of  the  computed 
impedances  using  the  simple  arc  and  the  final  form  versus  the  measured  data  is  seen  in 
Table  1  [6]. 

Figure  4  shows  the  calculated  current  distribution  on  the  LTWA  and  STWA  for  the 
22  kHz  test  case  (title  for  this  case  is  given  as  08:0322:90  VAL).  Both  amplitude  and  phase 
of  the  current  distribution  are  shown.  Labels  on  the  plot  list  the  various  aircraft  parameters 
used  in  the  calculation.  Analysis  of  the  results  for  the  amplitude  of  the  current  distribution 
shows  that  the  amplitude  is  peaking  at  60  Amps.  This  peak  occurs  at  an  arc  length  of  about 
11, (XX)  feet,  which  is  about  half  of  the  distance  from  the  drogue  to  the  feed  point. 

Table  1.  STWA  Geometry  Effects  on  Impedance 


FREQUENCY 

(KHZ) 

RUN  TITLE 

IMPEDANCE  1 

MEASURED 

NEC 

STWA  ARC 

NEC 

STWA  MODEL 

22 

08:0322:90V  AL 

600-j350 

712-j267 

647-j399.8 

WIRE  DC  RESISTANCE  =  4.5  OHMS/ 1000  FT 


Figure  5  compares  the  NEC  calculated  current  distribution  for  the  22  kHz  test  case 
with  a  pure  sinusoid^  current  distribution.  The  calculated  current  distribution  was 
normalized  to  a  maximum  value  of  1  for  comparison  purposes.  This  figure  demonstrates  that 
the  calculated  current  distribution  is  in  general  agreement  with  the  classical  sinusoidal  current 
distribution  approximation.  However,  there  are  some  slight  deviations  at  the  feed  point  and 
towards  the  drogue  end  of  the  LTWA. 

Table  2  is  a  comparison  of  measured  versus  calculated  input  impedance  for  the 
various  runs.  Figure  6  is  a  comparison  of  the  measured  and  calculated  input  impedances 
over  frequency.  Both  the  real  and  imaginary  parts  of  the  impedance  are  displayed.  The 
imaginary  part  is  negative  while  the  real  part  is  positive. 

The  calculated  impedance  data  agrees  well  with  the  measured  data.  The  close 
agreement  between  the  measured  and  calculated  impedances  provides  further  evidence  that 
the  NEC  model  is  closely  simulating  the  actual  dual  trailing  wire  antenna  system. 
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Table  2.  Measured  Versus  Calculated  Input  Impedance 


1  FREQUENCY 

RUN 

LTWA 

ltwa 

STWA 

MEASURED 

NEC 

1  (KHZ) 

TITLE 

DIST 

VERT 

VERT 

IMPEDANCE 

IMPEDANCE 

1 

08:0322:90VAL 

31.41 

69  81 

30.66 

60a-j350 

647-j399  8 

22 

08:0322:90V  AL 

85,32 

70  76 

30  66 

600-j350 

654  8-j396  9 

18 

08:03 182:90V  AL 

40.42 

72  89 

29  67 

45(>-j525 

642-j525 

17 

37.30 

68.41 

27.73 

48a-j525 

665.3-j491 

19 

09:2419:90V  AL 

88.18 

70.99 

29.66 

600-j490 

665.4-j459  4 

20 

09:2420:90V(B24) 

46.23 

68  18 

WM 

650^)490 

674.6-j437.4 

20 

09:2420:90V(B24) 

55.87 

69.94 

650-j490 

690.4-j432.2 

20 

09:2420:90V(E48)09:2 

19  11 

69.98 

30.28 

MSSM 

672.1-j433.1 

20 

420:90V(E48) 

73.81 

70  85 

30  28 

679  6-j430.6 

20 

09:2420:90V(E48) 

83.97 

68.26 

30.28 

657.9-j438.9 

21 

09:2421 :90V  AL 

74.02 

66.78 

30.09 

670-j480 

649-j423.3 

One  of  the  options  provided  in  the  NEC  code  is  the  ability  to  model  lossy  wire 
antennas.  NEC  uses  an  exact  expression  for  the  skin  effect  resistance  and  internal  reactance 
in  ohms/meter.  An  approximation  for  the  skin  effect  resistance  which  has  been  used  at 
moderately  low  frequencies  where  the  current  is  constant  throughout  the  wire  cross  section  is 
[7]: 


R  = 


1 

nroo 


(  \ 

4 

1 

''o 

^  48 

UJ 

(7) 


where 

ro  =  radius  of  wire 
a  =  wire  conductivity 
5  =  skin  depth 


At  HE  frequencies  a  better  approximation  can  be  given  for  both  the  real  and  imaginary  parts 
of  internal  impedance  of  a  round  wire  as: 


"HF 


(8) 


where 
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=  j_  = 


f  =  frequency 
/X  =  permeability 

A  more  accurate  expression  for  the  internal  impedance  of  a  round  wire  is  found  from 
evaluating  expressions  of  the  total  current  in  the  wire  and  the  magnetic  field  at  the  surface  of 
the  wire.  The  magnetic  field  can  be  obtained  from  the  electric  field  via  Maxwell’s 
Equations: 

V  X  I  = 


Expressions  can  be  derived  for  in  terms  of  Bessel  functions.  The  internal  impedance  per 
unit  length  can  then  be  found  from  the  ratio  of  the  tangential  component  of  the  electric  field 
at  the  surface  of  the  wire  divided  by  the  total  current  in  the  wire: 


Z,  =  - 


T/pCTro) 

2iirpOJo(Trp) 


where  T~ 


This  complex  internal  impedance  can  be  separated  into  real  and  imaginary  parts  using 
identities  of  Bessel  functions  resulting  in  Kelvin  functions. 

7  -  Ber  q  +  JBei  q 

‘  yj2nr^  [Ber'  q  +  jBei'  q 


where  q  = 


and  the  following  identities: 


Ber  v  +  jBei  v  =  Jq  ij'^^) 


(13a) 


Ber'  v  +  jBei'  v  =  —  {Ber  v  +  jBei  v)  =  j~'^Jl(j'^) 

dv 


(13b) 
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NEC  incorporates  these  exact  expressions  in  the  determination  of  wire  conductivity  skin 
effect. 


Table  3  is  a  comparison  of  NEC  calculated  input  impedance,  loss  resistance  and 
efficiency  for  various  antenna  configurations.  The  first  configuration  listed  is  for  a  half-wave 
vertical  dipole  at  an  altitude  of  20,500  feet  over  perfectly  conducting  ground.  It  is  offset  fed 
in  the  same  manner  corresponding  to  the  STWA  and  LlSvA  of  example  08:0322:90  VAL. 
The  parameters  are:  LTWA  =  19,500  feet,  STWA  =  2680  feet,  frequency  =  22kHz,  and 
wire  resistance  =  4.50/1000  ft.  This  case  would  correspond  to  an  antenna  with  a  verticality 
of  100%  and  offset  fed.  Shown  also  is  the  same  vertical  dipole  offset  fed  assuming  perfectly 
conducting  wire,  lossless  case.  The  loss  resistance  can  therefore  be  computed  from  the 
difference  of  the  real  part  of  the  impedance  for  the  two  cases  as  well  as  an  efficiency  factor 
given  by 


^lossless 

^lotsy 


(14) 


Also  shown  are  results  for  a  center  fed  vertical  halfwave  dipole  which  is  the  same  antenna  as 
the  offset  case.  The  results  are  for  over  perfectly  conducting  ground  and  in  free  space  for 
both  lossy  and  perfectly  conducting  wire.  The  last  result  is  the  model  for  the  dual  trailing 
wire  antenna  system.  The  efficiency  of  the  22  kHz  test  case  dual  trailing  wire  antenna  was 
found  to  be  47%,  assuming  a  known  wire  DC  resistance  of  4.5  ohms/10(X)  feet.  This 
suggests  that  nearly  half  of  the  power  input  to  the  antenna  is  being  dissipated  as  heat. 


IV.  Conclusions 


This  paper  has  shown  the  application  and  validity  of  using  the  Numerical 
Electromagnetics  Code  for  the  modeling  of  the  VLF  aircraft  dual  trailing  wire  antenna 
system.  A  description  was  given  of  the  approach  of  performing  mechanical  modeling  of  the 
steady-state  shape  geometry  of  the  dual  trailing  wire  antenna  composed  of  the  LTWA  and 
STWA.  Software  has  been  written  to  convert  these  wire  shapes  into  NEC  geometry  input 
data  involving  tapering  of  the  wire  segment  lengths  to  achieve  high  accuracy.  Figure  7 
summarizes  the  procedures  and  various  computer  codes  which  have  been  used  to  obtain 
results  in  this  work. 


Comparisons  have  been  made  to  actual  experimental  measurements  of  antenna 
feedpoint  impedance  with  excellent  agreement  to  computed  results.  Efficiencies  have  been 
calculated  from  various  modeled  results  using  exact  formulations  of  wire  conductivity  effects. 

These  results  are  important  in  that  the  current  distributions  as  computed  can  be  input 
to  an  additional  code,  TWIRENEC,  which  computes  VLF  propagation  in  the  earth- 
ionosphere  waveguide  [8].  Additionally,  similar  modeling  can  be  performed  with  dynamic 
mechanical  models  of  the  wires  to  obtain  the  variations  of  all  antenna  characteristics 
(impedance,  currents,  and  radiation  fields)  versus  time,  corresponding  to  the  mechanical 
motions  versus  time. 
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Table  3.  Efficiency  of  Wire  Antenna  (22  kHz) 


ANTENNA 

CONFIGURATION 

IMPEDANCE 

LOSS 

RESISTANCE 

EFF(%)  ! 

HALF-WAVE 
VERTICAL  DIPOLE 
OVER  GROUND 
(OFF-SET  FEED) 

1150-j85 

(LOSSY) 

724-j46 

(LOSSLESS) 

426 

63% 

HALF-WAVE 
VERTICAL  DIPOLE 
OVER  GROUND 
(CENTER  FED) 

157-l-j70.8 

(LOSSY) 

99.3-l-j50.5 

(LOSSLESS) 

57.7 

63% 

HALF-WAVE 
DIPOLE 
(FREE-SPACE) 
(CENTER  FED) 

133.7-l-j68.1 

(LOSSY) 

76.6-l-j47.2 

(LOSSLESS) 

57.1 

57% 

DUAL  TRAILING 
WIRE 

647-j399.8 

(LOSSY) 

303.8-j386.0 

(LOSSLESS) 

343.2 

47% 

WIRE  DC  RESISTANCE  =  4.5  OHMS/ 1000  FT 
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Test  Case. 


Long  TfoBtig  Who  Antenna  JLTWAJ 


Side  View  Showing  the  NEC  Wire  Segmentation  Geometry  for  the  22  kHz 
Test  Case. 


Figure  3.  Steady-State  Dual  Trailing  Wire  Antenna  Geometry  for  a  22  kHz  Test  Case. 
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Figure  4. 


Figure  5. 


LTWA  AND  STWA  CURRENT  DISTRIBUTION 
File:  08:0322:90VAL  Freq:  22  KHz 


Amplitude  Phase 


LTWA  DisL:  31.41  ft.  LTWA  Lenglti:  19500  ft.  Imp.  Tabb:  B22 

LTWA VerL:  69.81%  SIWALe#;  2680ft.  Meas.lmp.;  eOO'PO 

STWA  Vert;  30.66%  Allude;  20500  1  Catlm(L;  647  -|400  0 


NEC  Calculated  Current  Distribution  Resulting  From  the  22  kHz  Test  Case  of 
Figure  3. 


LTWA  AND  STWA  CURRENT  DISTRIBUTION 
File:  08:0322:90VAL  Freq:  22  KHz 


Normalized  Amplitude  vs  Sine  Curve 


LTWA  Disl;  31.411  LTWA  Lenglti;  19500  1  Iit^.T^;  B22 

LTWA  VerL;  69.81  %  STWAlengtii;  2680  1  Meas.  Imp.;  600  -  1350  0 

STWA  Vert:  30.66%  AM;  20500  1  CaLIn^.;  647  -j400  0 


Normalized  NEC  Calculated  Current  Amplitude  From  Figure  4  Compared  to  a 
Sinusoidal  Current  Distribution. 
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TACAMO  ANTENNA  IMPEDANCE 
Measured  vs  Calculated 


-o  -250 

<D 

Q. 

£  -500 


Imaginary 


X  Measured 
O  Calculated 


Frequency  (KHz) 

Figure  6.  Measured  Versus  NEC  Calculated  Input  Impedance  (Real  and  Imaginary  Parts) 
for  the  Validation  Runs  of  Table  2. 


WIRE4NEC 


STWANEC 


f  NEC  DATA  INPUT 
Wire  Segmenlation  Geometry  y 


NEC3D 


f  TWIRENEC  DATA  INPUT 
Wire  Segmentation  Geometry y 
Current  Distribution  / 


NEC  DATA  OUTPUT 
Antenna  Analysis 


TWIRENEC 


Electric  Field  Strength 
(uV/m) 


Current  Distribution 
Input  Impedance 
Near  Fields 
Power  Budget 
Efficiency 
Gain 


Figure  7.  Flowchart  showing  the  Interrelationship  Between  the  Program  Modules. 
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ABSntACT.  Different  surface  integral  equations  are 
presented  for  two-dimensional  composite  objects.  The  objects 
consist  of  impedance  bodies  partially  coated  with  dielectric 
materials.  In  all  of  the  formulations,  the  impedance 
boundary  condition  is  applied  on  the  impedance  surface  to 
reduce  the  matrix  size  in  the  numerical  solution.  The 
integral  equations  are  reduced  to  a  system  of  linear 
equations  via  the  point  matching  technique.  Application  of 
the  point  matching  technique  is  straight  forward  for  two 
dimensional  objects.  Because  of  this  surface  discontinuities 
can  be  treated  easily  without  the  problems  encountered  when 
using  triangle  basis  functions  as  a  result,  consideration  of 
two-dimensional  objects  gives  a  clear  picture  of  the  accuracy 
that  can  be  obtained  using  these  formulations.  Two  of  the 
formulations  discussed  herein  overcome  the  problem  of 
internal  resonance.  The  numerical  solutions  are  verified 
either  by  comparison  with  the  analytical  solutions  for 
cylindrical  objects  or  by  applying  self  consistency  tests  for 
targets  without  analytical  solutions. 


1  INTRODUCTION 

Receatly,  interest  has  been  renewed  in  using  the  impedance 
boundary  conditions  (IBC)  in  the  solution  of  electromagnetic 
scattering  problems.  Use  of  the  IBC  can  simplify  the 
solution  of  the  many  complex  electromagnetic  problems  for 
which  it  is  valid.  However,  composite  objects  of  complex 
structure,  in  terms  of  both  material  type  and  geometry,  are 
difficult  to  treat  using  the  IBC  because  the  IBC  is  often  not 
valid  on  all  object  surfaces.  Using  the  exact  boundary 
conditions  to  solve  these  problems  is  uneconomical  and 
requires  complicated  programming  when  compared  to  a 
me^od  incorporating  the  IBC. 

The  IBC  is  a  valid  approximation  under  certain  conditions 
[1-4].  The  problem  of  extending  the  IBC  for  use  on 
surfaces  v«diere  use  of  the  standard  IBC  would  usually  be 
considered  questionable  has  bera  investigated  to  some 
extent.  Generalized  impedance  boundary  conditions  were 
pressed  in  [5]  and  [6]  for  this  purpose.  Unfortunately,  use 
of  these  generalized  impedance  boundary  conditions  comes 
at  the  expense  of  considerable  analytical  complications  and 


requires  specialized  researchers  to  treat  each  new  geometry. 

An  alternative  to  using  generalized  impedance  boundary 
conditions  is  to  use  the  usual  IBC  only  on  surfues  where  it 
is  valid.  Since  IBC  is  a  localized  approximation,  it  can  be 
used  on  surfues  where  it  is  valid  and  the  exact  boundary 
conditions  can  be  used  on  the  rest  of  the  object  [7].  In  a 
practical  sense,  the  IBC  can  be  used  for  any  surfu^e  type  for 
which  the  surface  impedance  can  be  determined.  In  cases 
of  rapid  spatial  variation  of  the  surface  impedance, 
knowledge  of  the  derivative  of  the  surface  impedance  is  also 
required.  In  numerical  solutions  the  surface  impedance  must 
be  slowly  varying  to  allow  for  proper  piecewise 
approximation  on  the  surface  segmentation.  The  idea  of 
using  both  the  impedance  boundary  condition  and  the  exact 
boundary  conditions  in  a  single  formulation  leads  to  a 
technique  that  is  both  accurate  and  efficient. 

In  this  paper  such  a  technique  is  implemented  for  two- 
dimensional  (2D)  scatterers.  Here,  four  different  surface 
integral  formulations  are  implemented  for  two-dimmsional 
problems.  The  method  of  moments  using  the  point  matching 
technique  is  then  used  to  reduce  the  integral  equations  to 
matrix  equations.  The  separation  of  the  two  transverse 
polarizations  in  the  2D  problem  and  the  ability  to  easily 
implement  point  matching  in  the  numerical  solution  of  the 
2D  problem  (and  so  to  treat  surface  discontinuities)  will 
lead  to  a  more  complete  understanding  of  the  limitations 
imposed  on  this  technique  than  is  possible  with  three- 
dimensional  implementations. 


2  FORMULATION 

Consider  the  general  geometry  of  a  two-dimensional 
scatterer  consisting  of  an  impedance  body  that  is  partially 
coated  with  dielectric  as  illustrated  in  Fig.  1 .  The  impedance 
body  has  known  surface  impedance  and  the  dielectric 
coating  is  linear,  isotropic,  and  homogeneous.  For  this 
geometry,  there  are  three  distinct  regions:  V2  constituting 
the  impedance  body,  characterized  by  surface  inqsedance  1^; 
Vg,  the  exterior  of  the  scatterer,  characterized  by  the 
permittivity  and  permeability  of  the  free  space  (eg,  Pg);  and 


V,,  the  dielectric  region,  characterized  by  the  permittivity 
and  permeability  (c,,  /i,).  The  excitation  is  an 

Eo-E4e*.  H(rH<+R< 

Vo 


Origilul  Problem 

Fig.  1  Geometry  of  the  Original  Problem 


f  are  all  axially  directed  and  the  equivalent  magnetic 
currents  hf,  M* ,  and  ht  are  circumferratially  directed. 
Using  this  fact  the  equivalent  currents  can  be  written  in  the 
following  form: 


7’  =  7/  t  on  S" 

(1) 

=  M,"  i  on  5’ 

(2) 

where  q  represents  -1- ,  or  d,  t  is  the  unit  vector  in  the 

z  direction  and  t  is  the  unit  tangent.  The  unit  tangmt  is 
defined  by  the  equation 

1  ^  t  X  A.  (3) 


electromagnetic  plane  wave  of  incident  fields  E  and  tf.  The 
total  electric  and  magnetic  fields  in  region  V,  are  denoted  by 
E,  and  H„  respectively.  In  Vj  the  total  fields  are  not  of 
interest  and  have  therefore  been  assumed  to  be  zero  for 
convenience.  The  region  V,  is  bounded  by  E,  the  boundary 
surface  between  V,  and  Vj,  and  S',  the  boundary  surface 
between  V,  and  Vo.  The  region  Vo  is  bounded  by  S*,  the 
boundary  surface  between  Vo  and  Vj  ,  and  S'.  The  normal 
unit  vector  A  on  the  surface  S'  points  into  the  region  Vj 
and  out  of  the  region  V,.  On  the  surfaces  S*  and 
points  into  region  Vo  and  out  of  regions  V,  and  V,. 


In  the  zero  field  regions,  the  constitutive  parameters  are 
taken  to  be  the  same  as  in  the  non  zero  field  region,  so  that 
the  equivalent  currents  in  Figs.  2a  and  2b  radiate  into  an 
unbounded  homogeneous  medium. 

To  formulate  integral  equations,  the  impedance  boundary 
condition  is  applied  on  the  surfaces  S  and  S*  and  the 
continuity  of  the  tangential  components  of  the  electric  and 
magnetic  fields  on  S*  is  enforced.  These  boundary 
conditions  are  expressed  as 


The  equivalence  principle  is  applied  to  create  the  two 

V.  %  « 

ff,)  on  S' 

(4) 

auxiliary  problems  shown  in  Fig.  2  [7],  The  first  (Fig.  2a) 

is  the  exterior  equivalent  problem  that  is  electromagnetically 

equivalent  to  the  original  problem  in  region  Vg.  The  second 

n.  %  (<  ^ 

«o) 

(5) 

(Fig.  2b)  is  equivalent  to  the  original  problem  in  region  V,. 

In  the  TM  case  the  equivalent  electric  currents  f,  J* ,  and 

1 

II 

(6) 

^,1 

on  S* 

e4b'.  ii4h<  Vo 


Eitcrkr  Eqilvikioe 
(•) 


MraflcUi  Vo 


lK(riorB(|«N«leac« 

(b) 


Fig.2  The  equivalence  problems 


and 


A  X  H,  •  A  X  on  S"  (7) 

The  quantities  tiq,  q.  and  ri+  are  the  intrinsic  impedance  of 
the  free  space,  the  normalized  surface  inq>edance  on  S'  and 
the  normalized  surfKe  impedance  on  S*  respectively  (i;.  and 
ri+  are  normalized  by  ijg).  Equations  (4)  and  (S)  in^ly  that 
the  tangential  components  of  the  electrical  field  can  be 
expressed  in  terms  of  the  tangential  conqmnents  of  the 
magnetic  field.  In  terms  of  the  corresponding  equivalent 
magnetic  and  electric  surface  currents. 
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=  non.  71  X  <  =  non.-//  ? 


and 


•  non-  J'z  ^  A  =  %V.Jz  f 


PMCHW  formulation  [8]. 

The  currents  of  (1),  (2),  (8)  and  (9)  give  rise  to  TM  fields 
only.  Using  this  fact  equations  (10),  (11),  (6),  and  (7)  can 
be  reduced  to  scalar  integral  equations  and  written  as 


To  account  for  different  formulations,  the  integral  tions 
obtained  from  (4)  and  (S)  may  be  rewritten  as  follows  [7]: 

-  -  0  n-  (/»  X  «,)  onS-  (10) 

^0 


-  fiolu.  * n.  (0  X  «.)  O/IS*  (11) 
Vo 


where  a  and  0  are,  respectively,  the  combination  parameters 
weighing  the  EFIE  and  MFIE  just  inside  the  surfaces  S*  and 
5*.  Thus  different  field  formulations  can  be  obtained  by 
using  (10)  and  (11),  with  different  selections  of  a  and  0, 
together  with  (6)  and  (7).  Equations  (6)  and  (7)  represent 
the  PMCHW  boundary  condition  on  S‘.  These  formulations 
can  be  obtained  according  to  Table  I. 


-l£u  un*v.Vo  K 
Vo 

-0vAHM*v.Voffu  (A’f)  -  Hu  Ut)*H„ 


=  0  on  S' 


(12) 


--(fa.  U:)*V,Vo  Eoz  (Jz  i)*E^Uz")*Eoz(M.^] 
Vo 

*0vAH^U:)*V.VoHo.  Ut)*H^ 


^  5L  El-0n^H'  on  S' 

Vo 


(13) 


-lEo,  Uz^^V.VoEo.  (J:i)*E,^Ur>  ^V.vAz(J^'f> 

Vo 

*E„Uf)*E,^  Ut)  *  E^ 

=  —  E[  on  S* 

Vo 


Table  I 


Generation  of  different  formulations 


Formulation  tvne 

a 

& 

1. 

IBCE'PMCHW 

1. 

0. 

2. 

IBCH-PMCHW 

0. 

-l./rjc 

3. 

IBCC-PMCHW 

<  1. 

-I./tJc 

4. 

IBC-PMCHW 

1. 

1. 

In  Table  I,  stands  for  either  or  i\..  In  the  Erst  two 
formulations,  IBCE  and  IBCH  imply  that  E-  and  H-Eeld 
boundary  conditions  are  applied,  respectively,  just  inside  the 
impedance  surface  with  the  implementation  of  the  IBC 
approximation  for  the  magnetic  current,  whereas  PMCHW 
implies  that  the  continuity  of  the  tangential  Eeld  components 
is  enforced  on  the  dielectric  surface  In  the  third 
formulation  IBCC  denotes  the  combination  of  IBCE  and 
IBCH  on  the  impedance  surface.  In  the  fourth  formulation, 
IBC  implies  that  the  IBC  is  applied  on  the  impedance 
surface.  Indeed,  still  other  formulations  may  be  obtained 
when  the  Muller  formulation  is  applied  on  S'  instead  of  the 


Ho.  U:)*V,Vo  Ho,  (J:i)*HM*V.Vz/iuiJzi) 

*H„  Uz)^H^  (Mf)  -//„  <>5) 

=  Hi  on  S'' 

where 

H,  =  H  •  i  (1<) 

and 

E^  =  E  •  i  (17) 

The  operators  EJJ),  EJM),  HJJ),  and  HJ,M)  are 
determined  using  the  following  equations: 

E..  Ut)  =  -Jk,r,,  j^  j; ip')  g,(p,p')  dl  (18) 
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(Mf)  •  -2  •  V  X  W?  (pO  g,{p,  p')  dl>  (19) 


Ut)  -i‘  V  X  J\  (p' )  g,{p> 


HJMD-j^i-VxVx^^  Ml(p')g,(p,p')  dl'  (21) 

where 

g,(fi.t>')  -  1  H°'{k,\p-fi'\)  .  (22) 

V 

I  and  i},  =  .  In  the  expressions  above 

H,  and  6,  are  the  permeability  and  permittivity  of  region  V 

and  //o®  is  the  Hankel  function  of  second  type  and  zero 
order.  In  Equations  18  through  21  the  contour  integral  is 
evaluated  on  the  contour  that  results  from  the  intersection  of 
S*,  5^  or  with  the  x-y  plane.  These  contours  are  referred 
to  as  5^.  5*  or  5^  in  the  context  of  a  two-dimensional  body. 

The  contour  integrals  proceed  in  the  f  direction.  The 
vector  p  is  a  vector  in  the  x-y  plane  that  identifies  the  field 
point.  The  vector  p'  identifies  the  source  point. 

Equations  (12)  to  (IS)  are  specific  to  the  TM  incident  wave 
case.  In  the  TE  case  the  equivalent  electric  currents  T',  J* , 
and  /  are  all  circumferentially  directed  and  the  equivalent 
magnetic  currents  M‘,  M* ,  and  M'  are  axially  directed.  The 
TE  polarized  case  is  completely  dual  to  the  TM  case. 

The  scattered  field  can  be  calculated  in  the  exterior 
equivalent  situation  from  the  currents  J* ,  J*,  and  In 
equation  form  this  is  expressed  as 

£/(p.0) =n.  ’lo  (23) 

3  NUMERICAL  SOLUTION 

To  solve  the  surface  integral  equations,  the  contours  of  the 
scattering  body  are  divided  into  a  number  of  linear  zones. 
The  end  points  of  the  zones  lie  on  the  actual  contours  of  the 
body.  The  length  of  the  zones  is  taken  to  be  less  than  one 
tenth  of  a  wavelength.  The  currents  are  expanded  in  pulses 
basis  functions  multiplied  by  to-be-determined  coefficients. 
The  point  matching  technique  is  employed  to  reduce  the 
integral  equations  to  a  system  of  linear  equations  following 


the  procedure  given  in  [9].  Using  point  matching  and  pulses 
basis  functions  allows  for  accurate  representation  of  the 
currents  at  surface  of  the  discontinuities.  Accurate 
representation  of  the  currents  is  particularly  important  at  the 
junction  [10]  betw  een  the  dielectric  surface  and  conducting 
or  impedance  surface.  Equations  (18)  to  (22)  are  reduced 
to  a  standard  matrix  element  form  and  are  placed  in  the 
proper  location  in  equations  (12)  to  (IS)  to  obtain  the 
moment  method  matrix.  The  solution  of  the  matrix 
determines  the  current  coefficients  on  all  of  the  surfaces  of 
the  scattering  body.  These  coefficients  are  then  used  to 
obtain  the  far  scattered  fields.  The  expressions  used  to 
calculate  the  matrix  elements  and  the  scattered  fields  are 
given  in  [9]. 

4  RESULTS  AND  DISCUSSION 

The  numerical  solutions  of  the  four  formulations  defined  in 
Table  I  are  verified  in  this  section.  First,  circular 
cylindrical  bodies  are  considered.  The  series  solution  of  an 
impedance  cylinder  coated  with  a  linear  and  homogeneous 
dielectric  layer  of  uniform  thickness  is  used  to  verify  the 
numerical  results.  Fig.  3  shows  the  normalized  bistatic 
scattering  width  (ff/X)  from  an  impedance  cylinder  of  17. =0.5 
and  ka=3  coated  with  a  dielectric  layer  of  6, =4.0  and 
#*,=  1.0,  from  ka=3  to  kb=4.  The  agreement  between  the 
numerical  solution  of  all  the  formulations  and  the  series 
solution  is  satisfactory  for  both  TM  and  TE  polarizations. 
The  IBCE-PMCHW  and  IBCH-PMCHW  formulations  will 


<po 

Fig.  3Bistatic  scattering  width  of  a  coated  impedance 
circular  cylinder,  ka=3,  kb=4,  €,=4,  #*,=  1,  and  ij.=0.5. 

fail  when  the  impedance  core  is  at  resonance.  To  illustrate 
this  failure,  one  specific  example  is  presented.  The  case  of 
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a  circular  impedance  cylinder  coated  with  a  uniform 
dielectric  layer,  with  ka=2.71,  kb=3,  «,=2,  n,=  \,  and  ri. 


Fig.  4  Bistatic  scattering  width  of  a  coated  circular 
cylinder,  ka=2.71,  kb=3,  s,=2,  /i,=  l,  and  i}.=0.5. 

=0.5,  is  considered.  This  example  was  selected  because  it 
has  the  same  resonance  frequency  for  both  polarizations. 
Fig.  4  shows  the  conquirison  between  the  bistatic  patterns 
computed  numerically  using  the  present  method  and  the 
bistatic  pattern  determined  using  the  series  solution.  Clearly 
the  solution  due  to  the  IBCE-PMCHW  and  IBCH-PMCHW 
formulations  fails  to  give  the  correct  solution.  The  results 
from  the  other  two  formulations,  the  IBCC-PMCHW  and 
IBC-PMCHW,  are  in  good  agreement  with  the  exact 
solution,  indicating  that  they  are  not  affected  by  the  internal 
resonance  problem.  If  in  the  previous  example  the  surface 
impedance  is  zero,  the  perfect  conducting  core  case,  the 
IBC-PMCHW  formulation  will  fail  to  give  the  correct 
solution  because  it  reduces  to  the  E-PMCHW  formulation. 

In  the  previous  exanq>Ies  we  have  considered  completely 
coated  objects.  In  the  following  examples  we  will  consider 
partially  coated  objects.  The  bistatic  scattering  width  of  a 
half-dielectric/half-impedance  cylinder  of  ka  =  3,  i;.  = 

=  O.S,  e,=  l,  and  /t,  =  1  (’phantom”  dielectric)  is 
computed  for  all  formulations.  Results  for  the  object  with 
the  'phantom*  dielectric  half  must  be  the  same  as  those  for 
an  inqjedance  half-cylinder.  In  Fig.  S  a  comparison  is  made 
between  the  numerical  solution  of  the  impedance  body  with 
the  ’phantom”  dielectric  half  and  the  numerical  solution  of 
the  impedance  half-cylinder  with  i;,.  =  O.S.  Excellent 
agreement  is  observed.  The  same  geometry  is  considered 
with  6,=4  and  the  numerical  results  are  conqiared  with  the 


numerical  solution  of  the  same  object  where  the  impedance 
half  is  rq>iaced  by  a  lossy  dielectric  half  cylinder  (c,,  =  8-j  16 


Fig.  S  Comparison  between  the  bistatic  scattering  width  of 
a  half  impedance  cylinder  and  half  impedance/half 
dielectric  cylinder,  ka=3,  =0.5,«,=  1,  and 

and  ^,,=2-j4),  which  has  an  equivalent  normalized  surface 
impedance  of  0.5.  For  the  lossy  materials  the  exact 
boundary  conditions  are  enforced  on  all  the  object 
boundaries.  Fig.  6  and  7  show  a  comparison  between  the 
electric  and  magnetic  surface  currents,  respectively,  on  the 
outer  surfaces  of  the  object  using  the  numerical  solution 


MjPXW  phase  (deg.) 


S/  Xq  S/  Xq 
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Fig.  6  (TM)  Electric  outer  surface  current  on  a  half  lossy 
/half  dielectric  cylinder  and  half  impedance/half  dielectric 
cylinder,  ka=3,  e,=4,  and  ^,=  1 


IM,/E'I  phase  (deg.) 


0  0.5  1  1.5  2  2.5  3  0  0.5  1  1.5  2  2.5  3 


S/Ao  S/Xq 

Fig.  7  (TM)  Magnetic  surface  currents  on  a  half  lossy 
/half  dielectric  cylinder  and  half  impedance/half  dielectric 
cylinder,  ka=3,  e,=4,  and  /i,=  l 

with  exact  boundary  conditions  and  the  approxinute  solution 
of  the  E-PMCHW  formulation  for  the  TM  polarization. 
Good  agreement  between  both  results  is  obtained  in  the 
currents  magnitudes,  but  a  large  difference  is  observed  in 
the  phase  of  the  electric  current  on  the  pure  dielectric 
surface.  If  more  accurate  estimate  of  the  surface  impedance 
were  used,  better  accuracy  could  be  achieved.  A  significant 
observation  is  the  excellent  accuracy  of  the  surface  currents 
around  the  junction.  The  junction  between  three  or  more 
dielectric  regions  can  be  treated  accurately  even  when 
triangle  basis  functions  are  used  [10].  However,  the 
treatments  of  junctions  between  dielectric  and  conductor  or 
inqredance  surfaces  with  triangle  basis  fimctions  requires  an 
approximations  which  neglects  the  contributions  due  to  the 
magnetic  currents  on  half  the  triangle  basis  functions  around 
the  junctions.  This  problem  does  not  exist  when  point 
matching  and  pulses  basis  fimctions  are  used.  Fig.  8 
illustrated  that  the  exact  solution  and  the  approximate 
solutions  for  the  scattering  width  are  in  excellent  agreement 
with  each  other.  It  seems  that  the  current  error  has 
insignificant  effect  on  the  far  scattered  field  calculations.  If 
the  lossy  material  is  changed  so  that  =6-j4  and  /r,i  =  l, 
the  equivalent  surface  impedance  is  =0.3564 +j0. 108. 
For  this  example,  results  for  the  solution  incorporating  the 
exact  boundary  conditions  and  the  approximate  solutions  (the 
present  method)  are  compared  in  Fig.  9.  The  agreement  in 
this  case  is  not  as  good  as  in  the  previous  example.  This 
result  is  expected  because  the  surface  impedance  is 
calculated  assuming  that  there  is  no  wave  transmitted 
through  the  dielectric.  These  results  indicate  that  the  fields 
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Fig.  8  Bistatic  scattering  width  of  the  case  in  Fig.  6  and  7 


within  the  lossy  materia)  in  the  latter  example  are  much 
larger  than  in  the  former  example. 

Only  bistatic  radar  cross  section  has  been  considered  to  this 
point.  The  effect  of  different  angles  of  incident  can  be 
investigated  by  computing  the  monostatic  scattering  width. 
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Fig.  9  Bistatic  scattering  width  of  a  half  lossy/half 
dielectric  cylinder  and  half  impedance/half  dielectric 
cylinder,  ka=3,  «,=4,  and  /t,=  l. 
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Fig.  10  Monostotic  scattering  width  of  half 
impedance/half  dielectric  cylinder,  ka=3,  i}.=ij*=0.5, 
€,=4,  and 

The  monostotic  scattering  width  of  the  object  considered  in 
Fig.  8  is  given  in  Fig.  10.  The  agreement  between  the 
numerical  solution  of  the  different  formulations  is  within 
0.7S  dB  for  the  TM  polarization  and  0.25  dB  for  the  TE 
polarization. 

It  is  known  that  if  the  surface  impedance  of  an  object  equal 
to  the  intrinsic  impedance  of  free  space,  the  object  will 
have  zero  back  scattering  width.  This  suggest  that  a 
reduction  of  the  scattering  width  for  an  object  can  be 
obtained  by  manipulating  the  material  parameters  of  the 
object.  To  illustrate  this  technique,  the  dielectric-coated 
rounded  impedance  cylinder  shown  in  insert  of  Fig.  1  la  is 
considered.  For  the  original  scattering  object,  consider  the 
core  to  be  perfect  electric  conductor  (i;.=0.0)  coated  with  a 
uniform  dielectric  layer  of  thickness  t=0.  IXo  and  c,=4.0- 
jl.7.  If  the  transmission  line  nKxlel  is  used  to  calculate  the 
equivalent  surface  impedance  on  the  outer  dielectric  surface, 
the  surface  impedance  is  »j+  =  (0.6828+jl.0247)(neglecting 
the  curvature  of  the  surface).  To  reduce  the  scattering 
width,  the  core  is  selected  to  be  an  impedance  surface  and 
its  surface  impedance  value  is  manipulated  to  make  the 
equivalent  outer  surface  impedance  resistive  and  equal  to  the 
characteristic  impedance  of  free  space.  It  is  found,  using 
the  transmission  line  model,  that  with  % = (0.092  -t-j0.21S)  an 
outer  surface  impedance  of  i;,^=(1.0026-t-j0.0058)  results. 
The  scattering  width  which  is  calculated  for  the  original 


Fig.  11a  Monostotic  scattering  width  of  a  coated  rounded 
impedance  cylinder  when  »;.=0  and  »;.=0.092-fj0.215, 
«,=4-jl.7,  and  /*,=  1  (TM). 


object  and  for  the  reduced  scattering  width  body  is  shown  in 


Fig.  1  lb  Monostotic  scattering  width  of  a  coated  rounded 
impedance  cylinder  when  i;.=0  and  Tj.=0.092+j0.215, 
€,=4-jl.7,  and  /t,=  l  (TIE). 

Fig.  11.  In  Fig  11,  two  models  are  used  for  each  of  the 
bodies  to  calculate  the  scattering  width.  The  first  model  is 
the  two  surface  model.  The  two  surface  model  considers 
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both  the  inner  impedance  surface  and  the  outer  dielectric 
surface.  The  second  model  is  a  one  surface  model.  The  one 
surface  model  treats  the  outer  dielectric  surface  as  an 
inqiedance  surface  with  equivalent  surface  impedance  t)^. 
Results  of  the  one  surface  model  are  denoted  as  IBC  results 
in  Fig.  1 1 .  For  the  body  consisting  of  a  dielectric  coated 
core,  significant  reduction  of  the  scattering  width  is  achieved 
using  an  impedance  core  instead  of  a  perfect  electric 
conductor  (PEC)  core.  The  scattering  width  solution 
obtained  using  the  one  surface  model  is  less  accurate  than 
that  obtained  using  the  two  surface  model.  For  the  one 
surface  model,  the  solution  error  is  more  pronounce  in  the 
TE  case  than  in  the  TM  case.  This  example  illustrates  that 
proper  selection  of  the  core  material  results  in  low  back 
scattering  width  and  that  the  present  method  can  be  used  to 
accurately  predict  the  reduced  scattering  width. 

5  CONCLUSION 

Four  surface  integral  equation  formulations  are  developed 
for  two-dimensional  objects  composed  of  impedance  surfaces 
partially  coated  with  dielectric  material.  These  formulations 
are  useful  in  obtaining  accurate  and  economical  numerical 
solutions  for  bodies  that  can  be  modeled  as  coated 
impedance  surfaces.  Both  TE  and  TM  polarizations  are 
considered.  The  point  matching  technique  is  used  to  solve 
the  surface  integral  equations.  The  numerical  solutions  are 
verified  eiUier  by  comparison  with  the  series  solution  for 
circular  cylinders  or  by  comparison  with  exact  solutions 
(exact  boundary  conditions  on  all  surface  boundaries)  for 
other  objects.  The  internal  resonance  problem  is  investigated 
and  a  form  of  the  combined  field  integral  equation  is 
proposed  to  overcome  this  problem.  The  solution  accuracy 
is  shown  to  be  independent  of  the  polarization.  It  is  also 
shown  that  accurate  evaluation  of  the  surface  impedance  is 
very  important  to  achieve  good  accuracy.  One  example  is 
given  to  show  that  low  back-scattering  width  can  be 
achieved  by  selecting  the  material  properties  so  that  the 
outer  surface  impedance  is  equal  to  the  characteristic 
impedance  of  free  space. 
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Abstract 

The  spectral  domain  method  has  proved  to  be  a  suitable  analytical  tool  for  the  characterization  of  open, 
symmetrical  planar  structures.  The  method  requires  the  repeated  calculation  of  matrix  elements,  which 
each  involves  the  time-consuming  process  of  numerical  integration  over  an  infinite  range.  In  this  paper, 
suitable  basis  functions  for  the  expansion  of  electric  surface  currents  on  strips  or  electric  ftdds  in  slots, 
are  provided.  It  is  also  shown  how  the  use  of  these  basis  functions  makes  possible  the  efficient  and  rapid 
calculation  of  the  matrix  elements. 

1.  Introduction 

The  spectral  domain  method  has  become  an  important  tool  in  the  analysis  of  microwave  and  millimetre 
wave  integrated  circuits  [1].  Of  particular  interest  in  tiiis  paper  is  its  application  on  the  characterization 
of  open,  symmetrical  planar  structures  where  substrates  are  assumed  to  be  infinitely  wide  [2-6].  The 
method  generally  requires  a  significant  amount  of  analytical  preprocessing,  but  the  introduction  of  the 
inunitance  approach  [2]  has  simplified  the  derivation  of  the  dyadic  Green’s  function  elements.  However, 
some  practice  difficulties  are  still  encountered  during  the  implementation  of  the  method,  and  in  this  paper 
we  show  how  these  may  be  overcome. 

Application  of  the  immitance  approach  yields  the  spectral  dyadic  Green’s  function  for  the  planar  structure 
under  consideration.  The  unknown  electric  surface  currents  on  strips  (or  electric  fields  in  slots)  are 
expanded  into  finite  sets  of  basis  functions.  For  a  solution  of  the  dispersion  characteristics,  the  method 
requires  an  iterative  search  for  the  value  of  the  axial  propagation  constant,  0,  which  renders  the 
determinant  of  a  square  matrix  to  zero.  The  matrix  elements  n^  to  be  recalculated  during  each  iteration. 
The  computation  of  each  individual  matrix  element  requires  numerical  integration  over  an  infinite  range, 
where  the  integrand  contains  basis  functions  that  have  been  used  in  the  expansions.  These  calculations 
are  the  most  time-consuming  steps  in  the  implementation  of  the  method.  Due  to  the  slow  rate  of 
convergence  of  certain  integrals,  difficulties  are  encountered  in  attempts  to  attain  the  required  accuracy 
tolerances  during  the  numerical  integration.  In  this  paper,  we  provide  suitable  basis  functions  for  the 
expansion  of  unknown  electric  currents  or  fields.  The  use  of  these  basis  functions  facilitates  the  efficient 
calculation  of  the  matrix  elements. 

2.  Suitable  basis  functions 

Consider  the  general  multilayered  planar  structure  shown  in  Figure  1.  It  consists  of  a  number  of  dielectric 
layers,  with  strips  and/or  slots  spaced  symmetrically  about  the  y-axis  between  the  different  layers.  The 
structure  has  infinitely  thin  metallized  surfaces  witii  infinitely  wide  dielectric  substrates,  and  is 
homogeneous  in  the  z-direction.  It  may  be  analyzed  in  the  spectral  domain,  with  the  Fourier  transform 
defined  as 
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(1) 


♦  (a)  =  I  ♦  (jt)  ■'“*  dx 

•m 

«• 

$  (jt)  =  J_  f  ♦  (a)  « da 
2t  J 


The  electric  surface  currents  on  all  finite  width  strips  and  the  transverse  (with  respect  to  the  >  direction) 
electric  fields  in  finite  width  slots  are  expanded  into  finite  s^  of  basis  functions.  Using  the  inunitance 
approach  [2],  equivalent  circuits  may  be  constructed,  where  the  expanded  quantities  act  as  current  and 
voltage  sources  respectively.  From  Ae  equivalent  circuits,  the  elements  of  the  spectral  dyadic  Green’s 
function  are  obtain^  in  closed  form. 


X 


Figure  1  Symmetrical  multilayered  planar  structure. 


In  symmetrical  planar  structures,  the  strips  or  slots  in  a  specific  plane  y  =  yi  belong  to  one  of  two 
categories,  namely 

1 .  A  single  strip  or  slot  centred  at  x  =  0,  as  shown  in  Figure  2(a). 

2.  A  pair  of  strips  or  slots  centred  at  jc  =  b  and  Jt  =  -b  respectively,  as  depicted  in  Figure  2(b). 

A  multiconductor  transmission  line  supports  a  number  of  fundamental  modes.  For  symmetrical  structures, 
each  mode  belongs  to  one  of  two  mode  types,  which  we  denominate  as  even  and  odd  type  modes.  The 
field  distribution  of  even  (odd)  type  modes  is  characterized  by  the  fact  that  a  magnetic  conductor  (electric 
conductor)  may  be  introduced  in  the  plane  of  symmetry  without  disturbing  the  fields.  Unending  on  die 
mode  type,  only  basis  functions  which  are  either  even  or  odd  functions  of  jc  need  to  be  included  in  the 
expansion  of  currents  on  single  strips  or  fields  in  single  slots.  The  same  symmetry  considerations  hold 
for  the  case  of  a  pair  of  strips  or  slots.  It  is  however  necessary  to  include  basis  functions  which  are  both 
even  and  odd  with  respect  to  the  individual  axes  of  a  pair  of  strips  or  slots  [6].  This  could  be  explained 
by  means  of  the  following  example.  Let  Figure  2(b)  represent  a  pair  of  strips,  and  we  would  like  to 
expand  the  electric  surface  current  /,(z  ,y,)  for  an  even  type  mode.  From  the  symmetry  considerations 
it  then  follows  that  (x ,  y,)  should  be  an  odd  function  of  jc,  and  that  its  basis  (Unctions  should  thus  all 
be  odd  with  respect  to  jc  =  0.  However,  the  current  distribution  on  the  right  hand  strip  would  in  geno'al 
not  be  symmetrical  about  its  axis  ztx  =  b.  This  fact  therefore  requires  the  inclusion  of  basis  functions 
which  are  even  and  basis  functions  which  are  odd  with  respect  to  jc  =  6  in  the  domain  |jc  -  ^  w/2. 
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Figure  2  (a)  A  single  strip  or  slot  and  (b)  a  pair  of  strips  or  slots. 


The  surface  currents  or  electric  fields  in  the  plane  y  =  yj  may  be  expanded  using  any  kind  of  basis 
functions,  as  long  as  they  are  non-zero  only  for  |x|  ^  wl2  in  the  case  of  Figure  2(a),  or  only  for 
|x  ±  b  I  ^  w/2  in  Figure  2(b).  The  efficiency  and  accuracy  of  the  method  is,  however,  dependent  on  the 
choice  of  basis  functions.  The  singular  behaviour  of  the  electric  surface  currents  parallel  to  the  strip  edges 
or  the  electric  fields  normal  to  the  slot  edges  should  therefore  be  incorporated  in  the  basis  functions  [7]. 

In  general,  the  surface  currents  or  electric  fields  in  the  plane  y  =  y,  may  be  expanded  as 

=  E  aJAx)  C(x,y,)  =  E  ^  (2) 

n*!  «•! 


where  F  (j:  ,  y,)  and  G  (Jt ,  y,)  represent  /,  (jc ,  y,)  and  7,  (jc ,  y,),  or  (x ,  y,)  and  £,  (x ,  y,).  The  terms  a, 
and  b„  are  unknown  coefficients,  while /.(x)  and  g,(x)  are  basis  functions  which  satisfy  the  edge 
conditions.  In  the  spectral  domain,  this  becomes 

»*"1  «•! 


We  define  two  sets  of  functions  which  may  act  as  building  blocks  for  the  appropriate  bases.  These  are 
given  by 


7:..(2x/w) 
v/l  -  i2xlwy 


f„(x,w)  =^l-i2xfwy  f/„.,(2x/w) 


U|  S  w/2 
m  =  1 , 2,  3,. . 


(4) 


T„(x)  and  U„(x)  are  mth  order  Chebyshev  polynomials  of  the  first  and  second  kind  respectively.  These 
functions  are  shown  in  Figure  3  for  different  values  of  m.  {„(x ,  w)  is  singular  at  |x|  =  w/2.  Also  note 
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that  (t ,  w)  and  f„(jc ,  w)  are  even  functions  of  x  for  m  odd,  and  vice  versa. 


Their  Fourier  transforms  are  given  by 

Uo.w)  =7«- 


lja,w)  m  JJawl2)  (5) 


where  Jjix)  is  the  Bessel  function  of  the  first  kind  of  order  m.  The  basis  functions/,  (x)  and  g,  (x)  may 
thus  be  expressed  as  combinations  and  permutations  of  r«(x,H')  and  (.(x.w)  respectively.  Table  1 
provides  suitable  basis  functions  for  the  expansion  of  the  different  unknown  quantities. 


Single  /  Strip  /  Expanded 

Pair  Slot  quantity 


Single  Strip  y,(x,y,) 


Mode 

type 


Pair  Strip  J,(x,y,) 


Basis 

function 


{2«(JC,VV) 


^2n-,(X,W) 


f,(x+f>,»v)  - 
5  Ux-b,yv) 


Ux+b,w)  + 
5  Ux-b,w) 


ij^x-¥b,w)  + 

5  L^?^-b,w) 


^J^-hb,w)  - 
6 


S  L(x-b,yv) 


Ux+b,w)  + 
8  i„{x-b,w) 


5  Ux-b,w) 


8  Ux'b,w) 


Fourier  transform  of 
basis  function 


L(a,H') 


?a.-i(«,w) 


L-i(o.H') 


-2j  sin(ai>) 

2  cosiab)  iXfxM 


2  cos(af>)  f,(a,w) 
-2j  sin(of>)  f.Ca.w) 


2  cos(a&)  i,(a,H') 
-2/  sin(cKZ>)  |,(cK,w) 


-2y  sin(a^>)  |,(a,w) 
2  cos(a&)  L(a,w) 


-2y  sin(aZ>)  h(a,w) 
2  cos(ai>)  |,(a,w) 


2  cos(aft)  j,(a,w) 
-2y  sin(aZ>)  {,(a,w) 


2  cos(aft) 
-2/  sin(af>) 


-2/  sin(a6)  h(«.w) 
2  cos(aft)  l^ia,w) 


n  odd 
n  even 


n  odd 
n  even 


n  odd 
n  even 


n  odd 
n  even 


n  odd 
n  even 


n  odd 
n  even 


n  odd 
n  even 


n  odd 
n  even 


Table  1  The  nth  basis  function  for  the  expansion  of  currents  or  fields,  where  n=  1, 2, 3,  .  .  and 
6  =  ±  1  for  n  odd  and  n  even  respectively. 
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3.  Calculation  of  the  matrix  elements 


Application  of  Galerkin’s  method  together  with  Parseval’s  theorem  then  results  in  the  eigenvalue  equation 
which  needs  to  be  solved  numerically  [2-7],  The  matrix  elements  that  need  to  be  calculated  during  the 
process,  are  all  of  the  form  [7] 


P(j8)  =  I  G(a,0)  fi^(a)  da 

m 

=  2  I  G(a,0)  h^(a)  fi^(a)  da 
0 


(6) 


where  is  an  element  of  the  spectral  dyadic  Green’s  function,  hjia)  and  ht(a)  are  Fourier 

transformed  basis  functions  pertaining  to  a  domain  with  relevant  dimensions  denoted  by  subscripts  a  and 
b  respectively. 

For  closed  structures,  the  integrand  of  equation  (6)  might  have  poles  located  on  the  axis  of  integration. 
This  necessitates  the  application  of  residue  calculus  techniques  for  the  evaluation  of  the  integrals  [1,8]. 
These  processes  are  considered  to  be  beyond  the  scope  of  this  paper. 

In  the  case  of  open  structures,  no  poles  are  located  on  the  axis  of  integration,  and  therefore  no  special 
pole  extraction  techniques  are  required  for  the  calculation  of  the  integrals.  However,  the  matrix  elements 
are  computed  by  performing  numerical  integration  over  an  semi-infinite  range,  and  often  the  rate  of 
convergence  is  slow.  This  presents  difficulties  in  attaining  the  required  accuracy  tolerances.  By  applying 
the  following  techniques,  these  difficulties  may  be  overcome. 

Using  the  basis  functions  defined  in  (5)  and  Table  1,  we  see  that  the  product  hj[a)  htia)  consists  of : 

1.  A  complex  constant  term. 

2.  A  term  of  the  form  a  '  with  r  =  0,  I,  or  2.  For  example,  if  either  hjia)  is  proportionate  to 
h.(« .  or  is  proportionate  to  fja ,  Wj),  then  r  =  1 .  If  neither  hjia)  nor  hj(a)  is 
proportionate  to  then  r  =  0. 

3.  A  product  of  two  Bessel  functions. 

4.  Zero,  one  or  two  sin(aft)  and/or  cos(ab)  terms,  where  b  =  b„or  b  —  />*.  If  both  A,(a)  and  hjia) 
pertain  to  single  strips  or  slots,  no  sines  or  cosines  would  appear  in  the  product.  If  either  of  the 
basis  functions  pertains  to  a  pair  of  strips  or  slots,  then  the  product  of  the  bases  contains  one  sine 
or  cosine  term.  Finally,  if  both  hj^a)  and  ft»(a)  belongs  to  pairs  of  strips/slots,  each  would 
contribute  a  sin(a/i)  or  cos(a/i)  term  to  the  product  of  the  bases. 

When  the  product  of  two  trigonometric  terms  appear  in  the  integrand  of  equation  (6),  the  following 
identities  are  used  to  transform  it  to  the  sum  of  two  terms. 

sin(a/>J  sin(ot6j)  =  i{cosIa(i>„ -  cos[a(i>, +/>»)] } 

cos(a6J  cos(a*j)  =  i{cos[a(6, -/>*)]  + cos[a(ft^ +  (7) 

sin(a/>J  cos(a/>j)  =  i  {sin[a(/>,  - +  sin[a(ft^  ♦/>*)] } 


45 


Equation  (6)  may  thus  in  general  be  rewritten  in  the  form 


pm  =  E 

i*0 


m 

I 


J,(awJ2)  J,(awJ2) p,(a)  da 
of 


I 


=  E^. 

f.O 

/>o(«)  =  1 

/7,(a)  =  sin(obJ 

p,(o)  =  sin(ab*) 

p^(a)  =  cos(obj) 

p^(a)  =  cos[of(b^  +  b*)l 

^^(o)  =  sin(a(b,-b*)] 

P,(«) 

cos(oft,) 

sin[o(i», 

cos[ot(£>,-ft*)l 


(8) 


where  k  and  /  are  integers.  Ag,  i4,,  .  .  .<4,  are  complex  constants,  of  which  at  least  sevm  are  zero. 
Therefore,  a  maximum  of  two  terms  in  the  series  ne^  to  be  computed.  Note  that  equation  (8)  is  not  an 
expansion  of  the  original  integral  in  (6)  -  it  is  merely  a  general  representation  of  all  the  possible  forms 
equation  (6)  might  take  on  after  the  transformations  of  (7)  had  been  applied.  If,  for  example,  the  original 
integrand  in  (6)  contains  the  term  co&(abJ  cos(a6^,  the  only  non-zero  /4/s  in  (8)  are  and  A,.  For  die 
special  case  where  the  non-zero  constants  are  A^,  and  A^. 


For  t  ^  0,  the  integrands  of  the  F^03)  terms  are  all  oscillatory  over  the  entire  range.  These  integrals  have 
a  poor  convergence  rate  when  it  is  evaluated  with  conventional  quadrature  routines.  They  may  be 
computed  much  more  efficiently  by  treating  them  as  Fourier  integrals  with  the  sine  and  cosine  terms  as 
kernels.  A  special  routine  such  as  QDAWF  [9]  may  then  be  used  to  evaluate  these  integrals.  QDAWF 
is  an  adaptive  routine,  designed  to  integrate  functions  of  the  form  f(x)  sinfbuc)  or  /(x)  cosfcox)  over  a 
semi-infinite  range.  It  integrates  the  integrand  between  zeros  over  a  number  of  subintervals,  and  invokes 
an  extrapolation  scheme  in  order  to  estimate  the  integral. 

In  general,  the  Foifi)  term  may  be  calculated  by  using  a  routine  like  QDAGI 19].  QDAGI  is  an  integration 
routine  designed  to  numerically  evaluate  integrals  over  an  infinite  or  semi-infinite  range.  It  initially 
transforms  the  interval  into  the  finite  interval  [0, 1],  and  then  uses  a  21  point  Gauss-Kronrod  rule  to 
estimate  the  integral.  The  integrand  of  the  F„(fi)  term  is  well  behaved,  provided  that  the  basis  functions 
hjia)  and  h^ia)  are  not  associated  with  the  same  single  or  pair  of  strips  or  slots.  Inspection  reveals  that 
the  bases  and  h»(a)  then  necessarily  pertain  to  strips  or  slots  that  are  on  different  vertical  planes. 
(If  hjfii)  and  are  associated  with  different  strips  or  slots  on  the  same  vertical  plane,  the  constant  Ao 
is  zero,  and  therefore  the  FqOS)  term  need  not  be  computed.)  The  Green’s  function  elements  are  then  of 
such  a  form  that  the  integrand  decays  rapidly  for  large  values  of  a,  and  therefore  the  integral  converges 
quickly. 

However,  when  this  provision  does  not  apply,  the  basis  functions  hj(x)  and  hj(X)  are  defined  ovw  the 
same  domain  (so  that  w„  =  =  w  and  b„  =  b^  =  b).  The  integrand  then  decays  less  rapidly,  whidi 

causes  the  rate  of  convergence  to  be  slow.  As  reported  for  the  specific  case  of  a  single  slot  in  [4],  FqOS) 
may  then  be  converted  into  a  rapidly  convergent  integral  by  extracting  its  asymptotic  form  and  evduating 
it  in  closed  form,  lliis  is  done  by  stating  that 
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(9) 


F^m  =  I  /(a,/3)  da 
0 

•  m 

=  I  [/(a./S)-/.(«,/S)]rfa  *  \L(a,0)da 
0  0 

where  /.(a,/3)  is  the  asymptotic  form  of  the  integrand.  For  cases  where  convergence  is  slow,  the 
asymptote  is  usually  of  the  form 

=  COS)  — J,{awl2)J,(awl2)  (10) 


with  COS)  a  complex  function.  When  |it:  -  /|  =  2s  with  s  an  integer  (as  was  the  case  in  [4]),  the  second 
integral  in  the  final  expression  of  (9)  is  given  in  closed  form  by  [10] 


f 


f^(a,0)da 


cm  (-\y  I,(0wl2)  K,i$wl2)  k^l 
C(/S)(-iy/,(/Sw/2)^*(/3w/2)  k<l 


(11) 


where  I„(x)  and  KJx)  are  modified  Bessel  functions  of  the  first  and  second  kind  respectively. 


However,  when  (k  -  /(  =  2s  +  I  (i.e.  an  odd  integer),  this  integral  is  not  available  in  closed  form.  We 
therefore  need  to  modify  the  expression  for  /•(a,lS)  by  replacing  the  Bessel  functions  with  their 
respective  large  argument  forms  [11],  so  that 


/Jam  =  cm 


a 


- - - - —  COS(0(H') 

irow/2 


(12) 


and 


I  /Jamda  =  cm  (13) 


The  first  term  in  the  final  expression  of  equation  (9)  is  then  a  rapidly  convergent  integral,  and  is  suitable 
for  efficient  numerical  computation  with  a  routine  such  as  QDAGI  [9]. 

4.  Example 

Comparing  the  basis  functions  shown  in  Figure  3  with  the  conventional  trigonometric  bases  in  [3],  [6] 
and  [7],  shows  that  the  two  sets  of  functions  are  similar  in  form.  The  use  of  the  Chebyshev  bases  instead 
of  the  trigonometric  bases,  therefore  does  not  result  in  a  reduction  (or  an  increase)  in  the  number  of 
unknowns.  For  a  given  accuracy  level  of  the  final  result,  the  same  number  of  bases  should  be  included 
in  the  expansions,  irrespective  of  the  type  of  basis  functions  used.  When  the  special  calculation  techniques 
described  here  are  not  used,  the  CPU  time  required  to  compute  the  integrals  is  also  largely  insensitive 
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to  the  choice  of  basis  functions. 

However,  if  we  treat  oscillatory  integrals  as  Fourier  integrals,  and  apply  asymptotic  extraction  to  slowly 
converging  integrals,  the  required  CPU  time  is  reduced  significantly.  As  an  example,  we  consider  certain 
matrix  elements  that  are  calculated  during  the  analysis  of  the  semi  re-entrant  (SRE)  microstrip  section  16). 
This  structure  is  ideally  suited  to  illustrate  these  techniques,  since  it  comprises  of  both  a  single  strip  and 
a  pair  of  strips.  Instead  of  the  conventional  trigonometric  bases  used  by  the  authors  in  16),  we  now  utilize 
the  appropriate  basis  functions  specified  in  Table  1  to  calculate  the  following  elements  ; 

1 .  P 11  OS)  for  an  even  type  mode: 

From  the  definition  of  the  matrix  elements  (6,  eq.  (9)),  we  see  that 


Pn03)  =  I  2l!iaJ)l\(a)l\(a)da 

m 

=  2  I  2;; (a, ^ (a)  da 
0 


(14) 


where  2)i(a,/3)  is  an  element  of  the  spectral  dyadic  Green’s  function,  while  /I, (a)  is  a  Fourier 
transformed  basis  function  for  the  x  directed  current  on  the  single  strip.  From  Table  1,  it  follows 
that 7li(a)  =  f;(a  ,  w,),  which  yields 


I 

0 


a‘‘ 


J^{awJ2)  J^iawjl)  da 


(15) 


If  we  compare  this  to  equation  (8),  we  see  that  -  8  and  =  0  when  i  ^  0.  The  integral 
Fo(/S)  is  evaluated  as  indicated  in  equation  (9).  An  expression  for  COS)  is  obtained  by  calculating 
the  Green’s  function  element  for  a  >  >  /S,  so  that 


_ j 

+  /3^  «eo(* 


(16) 


and 


cm  = 


j 


(17) 


2.  Qll(fi)  for  an  even  type  mode: 

From  [6,  eq.  (9)],  we  see  that 

m 

Q,'2(/3)  =  f  2;/(a,^) />)/>)  (fa 

-OB 

g» 

=  2  J  2;/(a,^) />)/>)  da 

0 


(18) 
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Substituting  the  bases  yields 


Gj'IOS)  =  jllir^w^  I  ”  /,(aWj/2)  J^{awjl)  sin(oT)  da  (19) 


where  t  =  d  +  W2/2.  The  complex  constants  in  (8)  are  thus  given  by  ^,  =_/ 12  Wz  and  Aj  =  0 
when  { 1.  F,0?)  is  calculated  as  a  Fourier  integral. 

3.  5“  OS)  for  an  odd  type  mode: 

From  the  matrix  element  definition,  it  follows  that 


(20) 


which  reduces  to 


S”03)  =2ir^w^^  I  2^^(aJ)  J,{aw,l2)  J,iawJ2)  sinHar)  da 
0 
o» 

=  J  2”(a,/3)  J.iaw^n)  7,(aw,/2)  da 
0 

OD 

-  x^Wj^  I  2//(a,/3)  J^(aW2/2)  /j(aWj/2)  cos(2aT)  da 


0 

=  'Iq  ^o(i8)  *  ^6(13) 


(21) 


The  Ff,(fi)  term  is  treated  as  a  Fourier  integral  with  the  cosine  as  kernel,  while  the  integral  FqCS) 
is  calculated  by  extracting  its  asymptotic  form.  The  latter  is  done  by  noting  that  for  a  >  >  /S 


a 


"«o(^2*«r3) 


(22) 


so  that 


cm  = 


7/3^ 

««o(«r2"-«r3) 


2 


(23) 


These  matrix  elements  were  calculated  on  a  Persetel  PS8/90-3  computer,  in  the  one  case  by  utilizing  the 
special  techniques  described  in  section  3,  and  in  the  other  case  by  simply  integrating  the  integrands  as 
defined  in  equations  (1 4),  (1 8)  and  (20).  The  calculations  were  performed  for  an  SRE  structure  with 
dimensions  w/Xo  =  O.Ol,  Wj/)^  =  0.02,  dl\  =  0.005,  t/Xo  =  0.01,  s/Xq  =  0.0025,  e^  =  e^  =  2.2,  and 
with  V/3  =  X/Xo  =  0.7  .  The  integrals  were  calculated  with  a  relative  accuracy  criterion  of  0.(X)1  (i.e.  the 
error  should  be  smaller  than  0.1%  of  the  absolute  value  of  the  final  result).  The  CPU  times  required  to 
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calculate  the  matrix  elements  using  the  two  approaches  are  shown  in  Table  2. 


Element 

CPU  time  utilizing  special 
techniques  (seconds) 

CPU  time  without  special 
techniques  (seconds) 

P1103) 

0.37 

6.72 

0.20 

0.93 

5,^05) 

0.74 

12.15 

Table  2  CPU  time  required  to  calculate  the  different  matrix  elements. 


For  these  examples,  the  calculation  of  elements  with  oscillatory  integrands  as  Fourier  integrals,  and  the 
application  of  asymptotic  extraction  reduce  the  CPU  time  by  factors  of  about  4.S  and  18  respectively. 

5.  Conclusion 

Suitable  basis  functions  for  the  expansion  of  unknown  electric  currents  or  fields  as  required  by  the 
spectral  domain  method  applied  to  a  symmetrical  planar  transmission  line,  have  been  provided.  We  have 
shown  how  through  using  these  basis  hinctions,  the  required  CPU  time  for  the  matrix  element  calculations 
may  be  reduced  appreciably.  Elements  with  oscillatory  integrands  are  treated  as  Fourier-type  integrals, 
while  asymptotic  extraction  is  performed  to  enhance  the  convergence  rate  of  integrals  with  non-oscillatory 
integrands.  Application  of  the  latter  technique  has  been  limited  to  structures  with  single  strips  or  slots, 
but  with  the  additional  information  furnished  in  this  paper,  this  procedure  may  now  be  utilized  during 
the  analysis  of  any  symmetrical  multiconductor  transmission  line. 
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PARALLEL  IMPLEMENTATION  OF  THE  NUMERICAL 
ELECTROMAGNETICS  CODE 
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ABSTRACT 

NEC2  is  a  10  000  line,  public  domain,  FORTRAN  IV  program  for  electromagnetic  analysis.  This  program  has  been 
adapted  for  use  on  transputer  networks  of  various  dimension.  The  FORTRAN  code  has  been  modified  and  translated 
into  OCXZAM  where  necessary.  Parallel  algorithms  were  developed  and  implemented  for  each  NEC2  function  in  order 
to  optimise  efficiency.  Execution  efficiencies  in  excess  of  80%  were  attained. 


1  INTRODUCTION 

The  Numerical  Electromagnetics  Code  NEC2  (Burke, 
1981a,  1981b,  1981c)  was  originally  intended  to  run  on 
Mainframe  computers  but  has  recently  been  potted  to  run 
on  personal  computers  (PQ.  The  main  problem  with 
running  NEC2  on  a  PC  is  that  it  is  exceedingly  slow  and 
insuffi^nt  memory  limits  the  size  of  structures  that  can 
be  analyzed.  Stellenbosch  University  (Le  Roux,  1988) 
compil^  NEC  to  run  on  a  single  T800 INMOS  transputer 
and  showed  that  this  ran  much  faster  than  the  PC  version. 
NEC2  was  later  extended  (Nitch  and  Fourie,  1990)  to  run 
on  a  fixed  16  transputer  network  with  only  two  of  the  main 
algorithms  rewrittoi  to  execute  in  parallel  The  following 
di^dvantages  were  apparent: 

•  The  matrix  was  returned  in  full  to  the  host  transputer 
which  meant  that  problem  size  was  limited  by  the 
memory  on  the  host  transputer. 

•  Many  algorithms  were  still  sequential 

•  The  transputer  network  was  of  fixed  dimension 

This  paper  presents  solutions  to  these  problems.  It  is 
acknowledge  that  the  transputer  (TSOO)  is  qiute  an  old 
chip  and  is  relatively  slow.  The  parallel  sdgc^thms  dis¬ 
cussed  in  this  paper,  however,  may  also  be  applied  to 
contemporary  distributed-memory  parallel  processing 
machines. 


A. 
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where  is  the  interaction  between  segment  i  BOdJ  and 
is  a  function  of  the  wire  geometry 
/,  is  the  current  on  segment  i ,  which  is  obtained  by 
solving  the  equation. 

is  the  excitation  on  segment/ which  is  calculated 
from  the  specified  sources. 

Once  the  currents  on  the  segments  are  found,  other 
electromagnetic  characteristics  such  as  electromagnetic 
(EM)  fields  may  be  found. 

The  sequence  of  possible  events  carried  out  by  the 
program  is  as  follows; 

•  The  structure  geometry  is  read  from  a  file 

•  The  Z-matrix  is  filled  by  calculating  the  interaction 
between  segments.  This  requires  N*  operations  each 
of  which  involves  numeric^  integration. 


2  BACKGROUND  TO  NEC 

Structures  are  modelled  by  wire  segments  with  optkxis  to 
incliKle  sources,  loads  and  networks.  These  structures  may 
be  analyzed  in  various  environments  (free  space,  ground 
etc).  E^ntially  NEC  calculates  the  interaction  between 
the  N  wire  segments  making  up  the  structure  and  hence 
obtains  an  NxN  matrix.  An  excitation  vector  is  then 
calculated  as  a  function  of  the  sources.  Solving  this  matrix 
equation  nelds  the  currmts  on  each  segment  in  the 
structure.  Mathematically  this  may  be  expressed  as : 


•  The  Z-matrix  is  factorized  which  requires  simple 
(^rations. 

•  The  E-vector  is  calculated  from  the  sources. 

•  Solve  for  currents  which  requires  simple  oper¬ 
ations. 

•  The  efiea  of  /  networks  are  found  by  /  solve  oper¬ 
ations  on  the  Z-matrix  and  modifications  of  the 
E-vector. 
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•  A  single  EM  field  value  is  calculated  by  summing  the 
effect  of  the  structure  currents  at  that  point.  This 
involves  MN  operations  for  M  field  points. 

NEC  allows  structure  symmetry  to  be  exploited  in  order 
to  reduce  computational  effon  for  filling  and  factoring  the 
matrix.  This  involves  alteration  to  the  fill,  factor  and  solve 
of  the  matrix  which  will  be  discussed  in  more  detail. 

3  PARALLEL  ARCHITECTURE  AND 
COMMUNICATION 

The  transputer  (INMOS,  1989)  is  a  single  chip  micro¬ 
processor  with  4  communication  links  for  direct  com¬ 
munication  to  other  transi>uters.  The  general  architecture 
used  for  the  parallel  implementation  is  shown  in  Figure 
1.  It  should  be  noted  that  this  architecture  may  be  gen¬ 
eralised  for  processors  that  have  more  than  4  communi¬ 
cation  links. 


Each  transputer  has  knowledge  of  the  network  size  and  a 
number  identifying  its  position  in  the  network.  From  this 
information  it  is  possible  to  deduce  the  network  inter¬ 
connections  for  any  processor.  All  unconnected  links  are 
ignored  when  executing  the  algorithms  described  below. 

3.1  Broadcast  from  host  to  network. 

All  transputers  wait  on  link  0  for  a  message.  Upon 

receipt: 

*  the  main  transputer  on  a  board  whose  identity  is 
(proc.id  REM  4  =  0)  sends  the  message  on  links  1 
to  3.  (The  variable proc.uf  is  the  number  identifying 
the  processor  in  the  network,  whilst  the  XEM 
function  calculates  the  remainder  of  the  division  of 
proc.id  and  4.) 

•  otherwise  if  the  processor  is  not  the  main  transputer 
on  a  board,  then  send  out  to  next  board  or  cluster  if 
connected. 


CLUSTER  0 


3.2  Send  from  all  transputers  to  host 

All  transputers  send  their  own  message  out  on  link  0. 
Transputers  listen  on  links  1  to  3  for  messages  which 
will  be  routed  through  them  and  passed  on  through  link 
0. 

3  J  Broadcast  from  any  transputer  to 
network 

The  broadcasting  transputer  sends  its  message  on  all 
four  links.  Other  transputers  redirect  the  message  in  the 
following  way: 

*  If  the  broadcaster  is  on  the  same  board  as  the 
receiver  then  only  redirect  to  other  boards  or 
clusters 


Figure  1:  Hie  transputer  network. 

The  reasons  for  choosing  this  architecture  are  : 

•  the  path  fiom  any  proossor  to  any  other  processor  in 
the  network  is  minimized  when  compared  to  other 
networks  investigated  (various  meshes,  hypercubes 
etc.). 

•  the  algorithm  controlling  the  communication  is 
simple. 

•  extending  or  reducing  the  network  dimension  is  easy. 
This  extension  or  r^uction  may,  in  most  cases, 
achieved  by  adding  or  removing  single  processors 
from  the  network. 

Three  general  communication  strategies  were  required: 

•  Broadcast  fix)m  host  to  network. 

•  Send  from  all  transputers  to  host. 


*  If  the  broadcaster  is  not  on  the  same  board  as  the 
receiver  then  the  receiving  processor  redirects  the 
message  to  the  rest  of  the  processors  on  its  board 
and  to  connected  clusters  (not  other  connected 
boards). 

4  IMPLEMENTATION  OF  PARAL¬ 
LEL  ALGORITHMS 

4.1  Matrix  filling 

Inherently  the  matrix  filling  requires  each  transputer  to 
have  knowledge  of  the  structure  and  the  environment 
in  which  it  is  situated.  This  information  is  broadcast 
from  the  host  to  the  transputer  network.  This  enables 
each  transputer  to  calculate  any  matrix  element  inde¬ 
pendently  of  other  transputers.  It  is  important  to  decide 
on  which  part  of  the  matrix  each  processor  should  filL 
The  following  points  require  consideration  : 


•  Broadcast  from  any  transputer  to  the  rest  of  the 
network. 


Each  processor  should  calculate  approximately  the 
same  number  of  matrix  elements. 
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•  Each  processor  should  calculate  that  part  of  the 
matrix  which  it  requires  for  later  operations.  The 
reason  being  that  the  matrix  occupies  a  large  portion 
of  the  memory  and  hence  it  is  difficult  to  reshuffle 
efficiently.  In  the  previous  implementation  (Nitdi 
and  Fourie,  1990),  the  matrix  was  returned  to  the 
host  for  reordering.  This  method  obviously  did  not 
make  efficient  use  of  the  distributed  memory. 

The  matrix  was  fiUed  by  tows  because: 

•  both  column  and  tow  distribution  reduce  the 
communication  overhead  during  matrix  factoring 
when  compared  to  the  overhead  when  the  matrix  is 
divided  into  blocks. 

•  later  factoring  of  the  cyclic  block  matrix  (as  a  result 
of  structure  symmetry)  requires  that  it  be  in  row 
form  to  reduce  communication  overhead. 

Wrap  mapping  of  rows  was  employed  for  load  bal¬ 
ancing  in  the  factoring  algorithm  since  the  diagonal 
elements  of  the  Z-matrix  were  generally  the  largest. 

4.2  Matrix  factoring 

The  matrix  in  NEC  is  solved  using  Gaussian  elimination 
with  back  substitution.  The  paraUel  implementation  of 
this  algorithm  was  based  on  the  algorithm  presented  in 
a  paper  by  Giest  and  Romine  (1988)  and  is  shown  in 
Figure  2. 


FOR  k  =  OTOn-l 

determine  pivot  row  using  parallel  search 
update  permutation  vector 
IF  (I  own  pivot  row) 
broadcast  pivot  row 
ELSE 

receive  pivot  row 
ENDIF 

FOR  (all  rows  i>k  that  I  own) 

Ijjt  .'=  ttjj/au^ 

FOR  j  =  k+I  TO  n-1 

ENDFOR 

ENDFOR 

ENDFOR 


Figure  2:  The  parallel  algorithm  presented  by  Geist 
and  Romine. 

Where  lisa  temporary  vector  housing  the  pivot  col¬ 
umn. 

Implementing  this  algorithm  on  a  network  of  transputers 
requires  that,  for  good  load  balancing,  the  rows  of  the 
matrix  are  wrapped  onto  the  processors  (i.e  for  a 
network  of  16  processors  the  first  processor  should  have 
rows  1,  17,  33  etc  the  second  processor  should  have 
rows  2, 18,  34  and  so  on).  The  reason  for  the  wrapped 
row  mapping  is  that  once  a  processor  has  operated  on 
the  pivot  column  (column  k),  row  k  is  mt  used  in  any 
calculation  for  the  completion  of  the  algorithm.  The 


work  load  of  the  processor  is  therefore  reduced.  Thus 
to  ensure  that  the  processors  each  have  equal  loads 
throughout  the  execution  of  the  algorithm,  a  wrapped 
row  mapping  is  employed. 

The  communication  between  ptxxxssors  during  the 
computation  involves  only  the  broadcast  and  reception 
of  the  pivot  column. 


43  Matrix  solve 

Consider  the  lower  triangular  linear  system 
Lx  -ft 

where  L  is  a  lower  triangular  matrix  of  order  n 

b  is  the  right-hand  side  vector  of  dimension  n 
X  is  the  unknown  solution  vector 

The  serial  solution  of  this  system  may  be  represented 
by  the  code 


FOR  i  =  lTO  n 
FOR  j  =  lTO  i-1 

b,  =  fti-x/,^ 

ENDFOR 

X,  =  ft/Lfi 

ENDFOR 


Parallel  matrix  solve  algorithms  have  been  developed 
by  Guangye  and  Coleman  (1988),  and  Heath  and 
Romine  (1988),  and  others.  The  matrix  solve  routine 
used  in  this  implementation  was  based  on  an  algorithm 
presented  by  Heath  and  Romine  (1988)  and  is  shown 
below : 


FOR  j  =  lTOn 
IF  (I  have  row  j)  THEN 

eSdiP'^” 

fan-out  (x^  mcp(j)) 

FOR  (all  rows  i>j  diat  I  own) 

bi  =  bi-xjLii 

ENDFOR 

ENDFOR 


The  function  map(j)  relates  a  processor  to  the  tow  j. 
Thus  the  line  fan-out  (Xp  nu^>(j))  sends  the  message  Xj 
to  the  processor  with  row  j. 


4.4  Field  calculations 

Given  below  is  the  sequential  code  used  in  NEC  for  the 
far  field  calculation. 
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FOR  phi  =  1  TO  noOfPhiPoints 
FOR  theta  =  1  TO  noOfThetaPoints 
FORi  =  lTON 

calculate  field  at  (phi,tbeta)  due  to  current 
i 

vectorially  add  field  at  (phi,  theta) 
ENDFOR 

write  out  field  at  (phi,  theta) 

ENDFOR 

ENDFOR 


At  the  end  of  a  single  pass  through  the  inner  loop,  the 
results  are  written  to  disk. 

4.4.1  Field  calculation  code  for  the 
network  processors 

There  are  a  number  of  ways  one  can  split  up  the  nested 
FOR  loops  for  execution  in  parallel  Either  the  phi, 
the  theta  or  both  loops  may  be  divided  amongst  the 
processors.  The  number  of  theta  and  number  of  phi 
points  are  not  necessarily  the  same  for  each  radiation 
pattern  request.  Thus  splitting  the  phi  or  theta  loops 
could  produce  very  poor  perf^ance  figures. 

The  method  that  is  employed  finds  the  total  number 
of  radiation  points  for  the  field  calculation  and  hence 
reduces  it  to  a  single  loop.  These  points  are  divided 
amongst  the  processors  and  each  processor  finds  the 
field  at  these  points.  Implementing  this  method 
requires  the  decoding  of  radiation  pattern  point 
number  to  the  (theta,  j^i)  point  in  space. 

FOR  points  =  1  TO  noOfPhiPoints*noOf- 
ThetaFoints 
calculate  phi 
calculate  theta 
FORi=^lTON 

calculate  field  at  (phi,theta)  due  to  current 


I  vectorially  add  field  at  (phi,  theta)  | 

ENDFOR 

write  out  field  at  (phi,  theta) 

ENDFOR 


Splitting  the  loop  amongst  transputers  in  order  to 
achieve  parallel  execution  has  some  difficulties.  The 
transputers  in  the  network  do  not  have  access  to  disk 
and  memory  will  be  wasted  if  all  the  fields  are  stored. 
The  loop  must  hence  be  further  subdivided  such  that 
a  specified  number  of  fields  are  calculated  before 
infonnation  is  relayed  to  the  host. 

An  algorithm  employing  this  subdivision  was 
developed  for  the  transputer  network. 


4.4.2  Field  calculation  code  for  the 
host  processor 

The  principle  behind  the  algorithm  for  the  host 
processor  is  as  follows : 


FOR  itt  1  TO  noOfGroupsOfPoints 
request  network  to  find  radiation  pattern 
for 

ooiu)  of  points 

FOR  p^  ITO  noOfProcessors 
FOR  j  =  ITO  noOfPointsPerProcessor 
receive  radiation  pattern 
write  result  to  disk 
ENDFOR 
ENDFOR 

ENDFOR _ 

The  deficiency  in  this  algorithm  is  tb^  it  has  two  main 
serial  components.  The  first  is  the  request  to  the 
network  to  find  the  radiation  pattern  for  a  group  df 
points  and  the  second  is  the  writing  of  the  results  to 
disk.  Thus  the  network  waits  for  the  host  to  write  the 
results  to  disk  before  computing  the  next  set  of  fields. 

This  serial  emnponent  can  be  masked  by  buffering 
the  radiation  patterns  received  fitxn  tte  host.  A 
request  to  find  the  next  set  of  radiation  pattern  points 
may  be  made  before  writing  the  present  results  to  disk. 
Thus  at  the  expense  of  memory,  the  parallel  execution 
can  be  sped  up.  An  algorithm  employing  such  a 
bufiering  scheme  was  implemented  in  the  parallel 
NEC. 


4.5  Networks 

Networks  are  evaluated  in  NEC  through  the  use  of  a 
small  network  matrix  with  dimension  equal  to  the 
number  of  networks. 

The  following  steps  need  to  be  performed  for  the 
networks  in  the  structure:- 


Generate  the  RHS  vector  of  the  network 
equation  (Step  a) 

Fill  network  matrix  (Stq?  b) 

FOR  i  =  ITO  noOfNetworks 

Solve  the  Z-matrix  to  get  a  modification  vector 
(Step  c) 

use  modification  vector  to  adapt  network  matrix 
(Step  ^ 

ENDFOR 

Factor  the  network  matrix  (Stqp  e) 

Modify  RHS  of  network  matrix  equation  (Stq>  f) 
Solve  network  equations  for  voltage  across  ports  of 
those  networks  without  voltage  soures.  (Step  gf 
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The  oomputatioaally  time  ooasuming  portioa  of  this 
solution  is  filling  the  network  matrix  (steps  b,  c  and  d). 
Steps  a  and  b  are  carried  out  by  the  host  while  the 
traiisputer  network  is  filling  and  Sactori^  the  Z-matrix. 
Step  c  rrauires  the  use  of  the  factorised  ^matrix  to  find 
the  modification  vector.  While  c  is  performed  in  parallel 
on  the  network,  d  can  be  made  to  execute  ooncuireotly 
on  the  host. 


where 

A. 

Ai  - 

A. 

At 

A- 

. 

. 

A. 

Zj2  ... 

Arranging  the  code  in  this  manner  enables  the  host 
processor  and  the  network  to  work  in  paralkL  First,  the 
host  processor  fills  the  network  matrix  while  the  net¬ 
work  fills  and  factors  the  much  l^er  Z-matrix.  Then 
the  network  finds  the  modification  vector  (using  the 
parallel  algorithm  for  solving  discussed  before),  while 
the  host  processor  uses  a  previously  computed 
modification  vector  to  update  the  network  matrix,  using 
this  technique  the  time  required  to  fill  the  network 
matrix  is  approximately  equal  to  the  time  required  to 
find  the  modification  vectors. 

The  remaining  steps  are  comparatively  fast.  There  is 
little  point  in  factoring  the  network  matrix  on  the 
transputer  network  since  the  matrix  is  generally  of  small 
dimension  (30x30)  and  the  efficiency  of  the  network 
when  fiictorising  such  matrices  is  low.  The  solution  of 
a  matrix  of  this  dimension  is  not  very  time  consuming. 


4.6  Cyclic  Block  matrices 

Memory  and  Computation  time  is  saved  when  the 
structure  being  simulated  is  symmetric.  The  time 
required  to  fill  the  matrix  is  reduced  since  only  the 
interactions  between  those  segments  in  the  first  sym¬ 
metric  section  and  the  structure  are  calculated.  Hence 
the  filling  routine  is  simplified  to : 


FOR  i- 1  TO  noOfSegs 
FOR  j=  ITO  noOjSegsInSymSection 
find  interaction  between  segments  i  and  j  Le 
Zit 

ENDFOR 

ENDFOR 


The  resulting  matrix  is  structured  as  shown  below. 


[Aj  A2  ...  AjJ  X 


A 

- 

m 

• 

• 

■ 

. 

. 

A. 

and  s  is  the  number  of  segments  in  a  symmetric 
section. 

The  solution  of  this  system  of  equations  is  accomplished 
using  the  following  steps: 

•  Each  submatrix  is  combined  using  the  formula 

A,  - 

kml 

where  A,  is  the  i*^  submatrix. 

5^  are  factors  calculated  according  to  the 
type  of  symmetry. 

M  is  the  number  of  submatrices 

•  Each  submatrix  is  factored. 

•  The  excitation  vector  is  filled  in  the  normal  fashion. 

•  Each  submatrix  equation  is  solved. 

•  The  resulting  solutions  are  combined  using 

M 

4  ■  X  the  currents  on  the  segments, 

where  fj  is  the  solution  to  the  i*  submatrix  equation. 

Execution  of  these  steps  in  paraltel  may  be  considered 
in  two  sections,  namely,  the  filling  and  the  solution  of 
the  submatrices. 

The  filling  of  the  matrix  is  split  iq?  by  asking  each 
processor  in  the  network  to  fill  specific  rows.  It  is 
important  that  processors  fill  rows  of  the  Z-matrix,  since 
the  formula  used  to  combine  the  submatrices  opeike  on 
the  rows  of  the  submatrices.  Thus,  once  the  matrix  has 
been  filled,  each  processor  can  combine  the  elements 
of  its  portion  of  the  submatrices  without  having  to 
communicate  with  other  processors. 

Solving  the  submatrix  equations  is  accomplished  by 
sequentially  &ctoring  and  solving  eadi  of  the  submatrix 
equations  on  the  network.  The  resulting  solutions  are 
u 

then  combined  using/j  >  2  to  give  the  curraits  on 

km\ 

the  structure. 

5  PARALLEL  PERFORMANCE 

In  assessing  the  performance  of  an  algorithm  on  a  network 
of  transputers,  it  is  useful  to  compare  the  time  taken  to 
complete  the  task  on  the  network  to  the  time  taken  on  one 
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processor.  Thus  in  gauging  the  performance  of  the  paralk! 
NEC,  the  efficiency  and  speedup  of  the  computation  were 
calculated. 

Efficiency  and  speedup  are  defined  as  follows  :• 


.  Ametakentocompletetaskonone processor  , 

(timetakentocoHq>letetaskoHp  processors)  x  p\ 


iSpeed  Vp^ 


timet<Aentocompletetadc<moHe  processor 
dmetakentocompkutaskonpprocessors 


where  the  time  taken  to  complete  the  task  on  p  processors 
is  made  up  of  the  time  spmt  communicatmg  between 
processors  and  the  time  S(mt  doing  the  computation. 
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Calculating  tbeeffickocy  and  speedup  of  the  parallel  NEC 
requires  the  times  for  simulating  problems  on  both  a  single 
processor  and  on  the  network.  Since  foil  use  of  the 
distributed  memory  is  used  in  the  simulations,  a  single 
transputer  with  the  same  amount  of  memory  as  Sie 
network  should  be  used.  However,  a  processor  with  this 
amount  of  memory  was  not  available.  It  is  possible  to 
predict  the  time  that  it  would  take  for  a  single  processor 
to  do  a  sim  ulation.  Thus  the  efficiency  and  S{wedup  graphs 
use  some  predicted  values. 

Graph  1  shows  the  times  taken  to  simulate  structures  of 
various  electrical  size  on  processor  networks  consisting 
of  4, 8, 12,  and  16  processors. 
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Graph  1  :  Hw  times  taken  to  simulate  structures  of 
varying  electrical  size. 

Graph  2  gives  the  efficiency  of  the  transputer  networks 
and  graph  3  gives  the  speedup  of  the  simulations. 


Graph  2 :  Tlie  efficiency  of  the  transputer  network. 
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Graph  3 :  The  speedup  of  varioos  transputer  networks 
for  structures  of  varying  electrical  size. 

5.1  Radiation  Pattern  Performance 


Analysis  of  the  performance  of  the  transputer  network 
when  calculating  radiation  patterns  is  difficult  since 
there  are  many  factors  influencing  the  performance. 
These  factors  include : 

•  the  number  of  radiation  pattern  points  requested. 

•  the  number  of  segments  in  the  structure. 

•  the  number  of  processors  in  the  network. 

•  the  number  of  points  returned  to  the  host  at  a  time. 

•  the  speed  of  the  disk. 


When  computing  a  large  number  of  radiation  pattern 
points  for  a  structure  consisting  of  a  few  segments  there 
IS  a  bottle-neck  at  the  disk  since  the  netwo^  calculates 


the  points  faster  than  they  can  be  output  to  disk. 
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The  perfonnance  of  the  far  field  algorithms  are  illus¬ 
trated  in  graphs  1  and  S.  Graph  4shows  the  time  required 
to  find  a  number  of  radiated  fields  on  struaures  of 
varying  dimension  on  16  processors.  Graph  5  gives  the 
efficiency  of  the  process. 
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Grairfi  4  :  Tlw  time  required  to  And  a  radiation 
pattern  on  a  16  processor  network. 
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Graidi  5  :  Hw  efficiency  of  the  radiation  pattern 
caknilation  on  a  16  processor  network. 


6  CONCLUSION 

The  speedup  attained  by  the  transputer  network  indicate 
that  it  is  possible  to  significantly  reduce  the  execution  time 
of  NEC  by  distributing  the  program  onto  a  network  of 
processor  and  executing  the  code  in  paialkL  Full  use 
the  distributed  memory  was  made  by  careful  oonskteration 
of  the  operations  to  be  performed  on  the  largest  data 
structure  (the  interaction  matrix). 
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Abstract 

The  ability  of  the  MININEC3  code  to  accurately  model 
complex  radiating  structures  was  previously  unknown 
because  computer  compiler  constraints  limited  the  size 
of  analyzable  configuration.  With  the  conversion  of  the 
code  from  BASIC  to  FORTRAN  however  it  is  now 
feasible  to  use  it  to  simulate  such  structures.  This  paper 
presents  the  results  of  a  validation  exercise  aimed  at 
assessing  the  code’s  suitability  for  modelling 
configurations  which  involve  hundreds  of  segments. 

1  Introduction 

The  MININEC  moment  method  code  is 
sometimes  considered  erroneously  to  be  merely  a  more 
compact  version  of  the  widely  known  Numerical 
Electromagnetics  Code  (NEC).  There  are  though 
significant  differences  between  the  algorithms  used  in 
the  two  codes.  Of  particular  importance  are  the 
expansion  and  testing  functions  used  in  the 
representation  of  the  current  distribution.  MININEC 
uses  a  modified  Galerkin  approach  with  pulse  expansion 
and  testing  functions  whereas  NEC’s  employs  constant- 
sine-cosine  expansion  functions  and  delta  testing 
fimctions.  In  certain  moment  method  applications 
however  it  is  desirable  to  use  a  code  where  the  same 
expansion  and  testing  fimctions  are  employed  in  the 
representation  of  the  current  distribution.  One  particular 
application,  of  interest  to  the  authors,  is  the  calculation 
of  the  socalled  characteristic  modes  (Harrington  1975) 
of  a  radiating  structure.  This  is  facilitated  by  the 
solution  of  a  weighted  eigenvalue  equation  involving  its 
generalized  impedance  matrix  [Z].  If  a  Galerkin  or 
modified  Galerkin  moment  method  code  is  used  in  the 
initial  determination  of  [Z\  then  reciprocity  between 
segments  is  enforced  and  the  matrix  is  symmetrical. 
This  is  highly  advantageous  as  it  ensures  that  the 


characteristic  modal  currents  and  their  associated 
eigenvalues  are  real,  giving  them  a  clear  physical 
meaning  and  hence  allowing  greatm'  insight  into  a 
structure’s  radiation  characteristics. 

The  results  presented  in  diis  paper  relate  to  a 
research  program  ccmcemed  with  the  investigation  of 
the  characteristic  modes  of  complex  vehicular  structures 
(Murray  and  Austin  1993)  at  HF  and  VHF  frequencies. 
This  involved  the  use  of  a  moment  method  code  with 
the  ability  to  model  such  configurations  accurately.  The 
need  also  for  identical  expansion  and  testing  functions 
effectively  elimiruted  the  use  of  NEC.  Hence  the 
MININEC  code  was  chosen  as  the  basis  of  the  work.  A 
useful  feature  of  NEC  and  MININEC  is  the  similarity 
of  geometry  definition  which  means  that  the  wide  range 
of  available  NEC  pre-  and  post-processors  could  be 
employed  with  miniiiuil  modification.  The  one  major 
drawback  of  using  MININEC  was  that  it  had  not 
previously  been  used  to  model  complex  geometrical 
configurations  because  it  was  written  in  BASIC  and 
conqriler  constraints  imposed  severe  segmentation 
limits.  Conversion  of  the  code  to  FORTRAN  however, 
as  described  by  Miller  (1989a),  effectively  removed  this 
constraint  but  its  ability  to  nnodel  conqrlex  structures 
accurately  was  unknown.  This  paper  addresses  this 
issue  by  showing  the  results  of  a  MININEC3  validation 
exercise.  Three  approaches  were  adopted  using 
appropriate  aiudytical,  NEC-generated  and  available 
experimental  data. 


2  Development  of  the  MININEC3  Code 

The  Mini-Electromagnetics  Code  (MININEC) 
was  initially  developed  by  Julian  er  al  (1982)  as  a 
Method  of  Moments  code  for  the  analysis  of  sinqrle, 
relatively  small  (in  terms  of  wavelengths),  anterma 
structures.  The  original  purpose  of  the  code  was  to 
provide  a  tool  for  the  rapid  analysis  of  such  sittqrle 
anteimas  on  small  micro  or  desktop  computers.  It 
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allowed  the  computation  of  a  number  of  the  basic 
characteristics  such  as  current  distribution,  input 
impedance  and  far-field  pattern  of  an  antenna  situated  in 
either  free  space  or  over  an  infinite,  perfectly 
conducting  ground-plane.  There  were  many  limitations 
to  the  initial  version  of  MININEC.  Point  reactive 
loading  was  available,  although  this  was  not  allowed  for 
wire  segments  intersecting  the  perfectly  conducting 
ground  plane.  Also  any  wire  that  did  intersect  the 
ground  plane  had  to  do  so  at  an  angle  of  90°.  In 
addition  computer  technology  of  the  time  introduced 
further  constraints.  Computer  memory  and  the  available 
BASIC  compilers  limited  the  size  of  analyzable 
geometries  to  around  30  segments.  This  meant  that  the 
maximum  size  of  antenna  that  could  be  handled 
confidently  was  of  the  order  of  1  wavelength. 

A  second  version  of  the  code  was  developed  by 
Li  et  al  (1983)  which  addressed  some  of  the  ground 
plane  constraints  listed  above.  The  next  major  step  in 
the  progression  of  the  code  however  was  the 
development  of  MININEC3  by  Logan  and  Rockway 
(1986),  together  with  the  rapid  development  of 
computer  technology.  Although  the  program  was  still 
written  in  the  BASIC  language,  compilers  were 
available  that  could  address  up  to  64k  of  memory, 
enabling  antenna  structures  of  up  to  12S  segmrats  or  of 
the  order  of  8  wavelengths  to  be  analyased.  The  relative 
increase  in  desktop  computer  speed  also  made  the 
analysis  of  such  structures  feasible  within  a  reasonable 
period  of  time. 

This  third  version  of  MININEC  also  included 
many  more  of  the  features  that  are  available  with  the 
mainframe  NEC  code.  As  with  the  previous  codes, 
current  distribution  and  hence  input  impedance,  along 
with  far  field  patterns  could  be  calculated.  The  feature 
of  point  loading  the  anteima  was  improved  to  include 
the  ability  to  add  either  fixed  or  frequency  dependant  s- 
domain  loading.  The  anteima  could  now  be  situated  in 
either  free  space  or  over  a  perfect  or  finite  conducting 
ground  plane.  This  feature  used  the  Fresnel  reflection 
coefficient  approximation  to  allow  five  changes  in 
ground  impedance  with  distance  using  either  circular  or 
linear  boundaries.  Near  fields  were  also  computed,  if 
required. 

Other  variants  of  the  MININEC  code  have 
bera  developed  that  have  improved  the  interface  with 
the  user  and  also  presented  improved  routines  for  the 
display  of  calculated  data.  Such  codes  are  described  by 
Lewallen  (1991)  and  include  The  MININEC  System 
(Logan  1988),  MN4  (Beezley  1992)  and  ELNEC 
(Lewallen  1991).  Although  not  improving  the 
capabilities  of  the  basic  algorithm  or  speed,  they  present 
an  attractive,  "user  friendly*  interface  and  rapid 
methods  for  displaying  calculated  data  graphically. 


The  algorithms  used  in  the  MIN1NEC3  code 
and  its  variants  provide  a  reliable,  stable  method  for  the 
calculation  of  generalized  anteima  characteristics.  There 
is  no  theoretical  limit  to  the  size  of  geometry  that  can 
be  analyzed  using  these  algorithms,  although  it  is 
limited  by  two  external  factors.  Firstly,  as  discussed 
above,  the  size  of  core  memory  that  a  BASIC  conqiiler 
can  address  limits  the  physical  size  of  the  problem. 
Secondly,  the  speed  of  operation  of  the  computer  on 
which  the  software  is  mounted  places  a  further  practical 
constraint  on  the  time  available. 

The  first  of  these  constraints  was  addressed 
with  the  implementation  of  a  version  of  MININEC3  in 
FORTRAN  as  described  by  Miller  (1989a).  The 
FORTRAN  programming  language  overcomes  all  the 
inherent  constraints  of  BASIC  compilers.  There  is  no 
theoretical  limit  to  the  size  of  memory  that  can  be 
addressed  and  the  greater  efficiency  of  the  language  for 
numerical  calculations  improves  execution  time.  Also 
the  code  may  be  used  with  a  wider  range  of  conqmters 
including  mainframes. 

3  Validation  of  MIN  IN EC3 

Three  differrat  approaches  will  be  used  in  the 
assessment  of  the  ability  of  M1NINEC3  to  model 
complex  antenna  structures.  Firstly  the  analytical  results 
for  the  input  impedance  of  a  V-dipole  antenna  are 
examined.  This  configuration,  with  an  acute  angle 
between  two  conductors,  has  been  reported  to  produce 
erroneous  impedance  results  with  some  moment  n^thod 
codes  (Austin  1993).  Since  wire  grid  models  of  conqilex 
structures  frequently  contain  such  acute  angles  it  is 
therefore  inqwrtant  that  the  ability  of  MININEC3  to 
model  the  situation  is  considered.  Secondly  MININEC3 
was  used  to  calculate  the  radiation  characteristics  of  a 
number  of  conqilex  wire  grid  models  representing 
continuous  surfaces.  These  are  compared  to  the  results 
obtained  using  NEC.  Finally  the  experimentally- 
determined  input  impedance  of  a  monopole  mounted  on 
a  conducting  box  is  conqiared  to  the  results  obtained 
with  MININEC3. 


4  The  V-dipole  Antenna 

The  input  inqiedance  of  a  centre-fed  V-dipole 
antenna  in  free  space  was  determined  analytically  by 
Jones  (1976).  His  results  showed  excellent  agreement 
with  experimentally  determined  data  for  a  base  fed 
monopole  antenna  positioned  at  various  angles  over  a 
large  conducting  ground  plane.  The  antenna 
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configuration  used  by  Jones  is  shown  in  figure  1.  He 
varied  the  angle  i'  between  the  two  wires  from  30*  to 
180*  in  30*  incremoits. 


For  the  MININEC3  model  three  uniform  and 
one  bqjered  segmentation  scheme  were  employed  and 
are  shown  in  table  I.  Schemes  1,  2  and  3  are  the 
uniform  segmentation  configurations.  Notice  that 
schemes  2  and  3  involve  higher  densities  than  the 
generally  accepted  MININEC3  "rule  of  thumb*  of  20 
segments  per  wavelength  of  wire.  Scheme  4  is  the 
graded  segmentation  case.  This  consists  of  dividing  each 
wire  into  two  equal  sections.  For  the  section  including 
the  feed-point  a  higher  segmentation  density  was 
employed.  Figures  2  and  3  then  show  the  predicted 
input  impedances  obtained  using  these  four  schemes  and 
compares  them  to  those  of  Jones. 


Table  I.  V-dipole  segmentation  schemes. 


Scheme 

1 

2 

3 

4 

Total  segments 

12 

20 

40 

30 

Segmentation  density 
(per  wavelength) 

24 

40 

80 

80- 

Region  1 
40- 

Region  2 

Figure  2.  Comparison  of  analytical  and 
MIN^1EC3  predicted  input  resistance. 


Figure  3.  Comparison  of  analytical  and 
MIN]^C3  predicted  input  reactance. 

Figure  2  shows  that  the  MININEC3  input 
resistances  obtained  with  all  of  the  segmentation 
schemes  are  virtually  indistinguishable  from  those  of 
Jones.  For  the  input  reactance,  shown  in  figure  3, 
however  it  is  noticeable  that  as  'F  is  reduced  the 
difference  between  MININEC3  and  the  analytical 
results  is  substantial  for  12  and  20  segments.  By 
contrast  schemes  3  and  4  are  in  excellent  agreement  for 
all  'F. 

It  is  apparent  from  these  results  that 
MININEC3  is  capable  of  accurately  modelling  wires 
with  an  acute  angle  provided  that  an  appropriate 
segmentation  scheme  is  employed.  The  results  here 
suggest  however  that  only  in  the  local  region  around  the 
junction  between  the  wires  is  it  necessary  to  use  this 
relatively  dense  segmentation  scheme.  Thus  for 
♦  <90°,  schemes  3  and  4  show  good  agreement 
compared  to  the  poor  agreement  of  schemes  1  and  2. 
The  advantage  of  scheme  4  compared  to  scheme  3  is  the 
reduction  in  computation  time.  Tapering  the  segments 
in  the  vicinity  of  the  feed  has  also  been  shown  by 
Lewallen  (1991)  to  improve  the  accuracy  significantly. 


5  Comparison  with  NEC 

The  NEC  program  is  the  one  of  the  most 
widely  used  and  highly  regarded  momeat  method  codes 
for  the  analysis  of  complex  antenmi  systems.  Comparing 
its  results  to  those  of  MININEC3  is  therefore  a  useful 
validation  exercise.  A  wide  range  of  continuous  surfaces 
represented  by  wire  grid  models  were  therefore 
modelled  using  both  codes.  This  section  compares  both 
the  predicted  far-tieid  patterns  and  input  impedances  of 
two  specific  geometrical  configurations.  The  first  is  a 
0.4X  X  0.6X  flat  plate  with  a  monopole  mounted  in  the 
centre  as  shown  in  figure  4.  The  second  is  a  more 
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complex  box  structure  containing  a  aperture  and  two 
vertical  monopoles  as  shown  in  figure  3.  Wire  grid 
models  of  these  structures  were  developed  using  the 
mesh  generator  of  Najm  (1991).  The  grid  spacing  is 
O.IX  and  the  so-called  equal-area  rule  was  applied  to 
determine  the  wire  radius.  Figures  6  and  7  show  the 
resulting  wire  meshes. 


z 


Figure  4.  Flat  conducting  plate  with  centre 
mounted  monopole. 


t 


Figure  5.  Conducting  box  with  aperture  and  two 
monopoles. 


Figure  6.  Wire  grid  model  of  the  geometrical 
conflguration  shown  in  figure  4. 


Figure  7.  Wire  grid  model  of  the  geometrical 
configuration  shown  in  fqjure  5. 

For  the  first  model  the  mcmopole  was  excited 
with  a  single  voltage  source  situated  at  its  base.  The 
two  monopoles  of  the  second  configuration  however 
allowed  them  to  be  base  fed  with  voltage  sources  of 
different  magnitude  and  phase.  This  is  clearly  useful  in 
a  validation  exercise  as  numerous  dissimilar  far-fields 
are  obtaiiuible,  as  well  as  complex  current  flow  on  the 
whole  structure. 

Considering  firstly  the  predicted  radiation 
patterns  in  free  space  of  the  simpler  plate  configuration, 
figure  8  shows  the  variation  in  elevation  of  the  vertical 
£-field  component,  while  figure  9  shows  the  azimuthal 
variation  of  the  horizontal  £-field  component.  Excellent 
agreement  is  achieved  in  each  case  with  a  maximum 
difference  of  less  than  IdB  at  any  point.  Next 
considering  the  more  complex  box  structure  figures  10 
and  1 1  show  the  azimuthal  variation  of  the  vertical  E- 
field  conqwnent  with  the  monopoles  fed  firstly  in-f^use 
and  secondly  in  anti-phase.  Figures  12  and  13  ^ow  the 
variation  in  elevation  of  the  horizontal  £-field 
component  for  the  same  feed  configurations.  Overall  the 
agreement  shown  for  this  configuration  is  reasonably 
good  with  the  patterns,  in  general,  exhibiting  the  same 
shape  with  nearly  idratically  positioned  peaks  and  nulls. 
The  most  noticeable  difference  in  the  pattern  plots  is  in 
the  depth  of  the  nulls  udiich  in  the  worst  case  is  up  to 
SdB.  Overall  though  both  codes  are  predicting 
substantially  similar  spatial  distributions  of  power. 

Next  we  consider  the  input  impedance  where 
table  II  shows  the  values  predicted  by  both  codes.  It  is 
noticeable  that  the  largest  differrace  between  the  two 
sets  of  results  is  in  the  reactive  componoit.  MININEC3 
consistratly  predicts  an  input  reactance  that  is 
substantially  more  capacitive  than  that  of  NEC.  Miller 
(1989b)  identified  a  similar  effect  >%diere  differmt 
computer  models  occasionally  produce  similar  results 
but  with  a  noticeable  frequency  shift  from  one  another. 
The  overall  dimoisions  of  the  structures  considered  here 
are  around  the  first  natural  resonance  where  generally 
the  reactive  componoit  of  the  input  impedance  is 
changing  rapidly. 
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(0=O“) 

F^ure  8.  Variation  in  elevation  of  the  MININEC3- 
and  NEC-predicted  vertical  £-fleld 
component  for  the  flat  plate. 


0  <0  120  JIO  240  300  300 

fkl  (dtfrttt) 
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Figure  9.  Azimuthal  variation  of  the  MIN1NEC3- 
and  NEC-predicted  horizontal  £-field 
component  for  the  flat  plate. 


(0=90”) 

Figure  10.  Azimuthal  variation  of  the  MININEC3- 
and  NEC-predicted  vertical  £-fieid 
component  for  the  box  with  aperture 
with  the  monopoles  fed  in  phase. 


(«=90“) 

Figure  11.  Azimuthal  variation  of  the  MININEC3- 
and  NEC-predicted  vertical  £-fleid 
component  for  the  box  with  aperture 
with  the  monopoles  fed  in  anti-phase. 


(<^=90“) 

Figure  12.  Variation  in  elevation  of  the  MININEC3- 
and  NEC-predicted  horizontal  £-field 
component  for  the  box  with  aperture 
with  the  monopoles  fed  in  phase. 


(^=90») 

Figure  13.  Variation  in  elevation  of  the  MININEC3- 
and  NEC-predicted  horizontal  £-field 
component  for  the  box  with  aperture 
with  the  monopoles  fed  in  anti-phase. 
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Table  II.  NEC  and  MININEC3  predicted  input  impedances. 


Geometrical  Configuration 

NEC 

MININEC3  1 

Flat  plate 

43.34+j42.43 

37.66 -(-jl8.08 

Box  with  cavity 
(fed  in-phase) 

Monopole  1 
Monopole  2 

151.64+j25.55 

57.10+j22.26 

148.88+j7.14 

39.80-jl.73 

Box  with  cavity 
(fed  in  anti-phase) 

Monopole  1 
Monopole  2 

66.46 +j21. 87 
36.59+jl5.21 

54.64+jl.28 

27.45-j4.98  | 

Table  ID.  NEC  and  MININEC3  predicted  input  impedances  with  an  empirically  determined  frequency  shift  in 
the  NEC  results. 


Geometrical  Configuration 

NEC 

Frequency  =  0.95^ 

M1NINEC3 

Frequency  =  ^ 

Flat  plate 

34.46 +j23.21 

37.66 -l-j  18.08 

Box  with  cavity 
(fed  in-phase) 

Monopole  1 

150.70+j9.21 

148.88-l-j7.14 

Monopole  2 

52.15+j7.26 

39.80-jl.73 

Box  with  cavity 
(fed  in  anti-phase) 

Monopole  1 

64.39-jl.42 

54.64-f-jI.28 

Monopole  2 

35.12-jl.34 

27.45-j4.98 

Hence  reducing  the  frequency  of  the  NEC 
models  by  an  empirically  determined  value  of  5  %  gives 
the  input  impedances  shown  in  table  III.  The  original 
MININEC  results  are  also  shown  for  comparison. 

The  frequency  shift  of  5  %  clearly  reduces  the 
difference  between  the  two  sets  of  inq)edance  results. 
Whereas  the  resistive  component  of  the  NEC  results 
generally  changes  negligibly,  the  NEC-predicted  input 
reactance  of  each  feed-point  is  reduced  considerably. 
Good  agreement  is  now  achieved  between  the  two  sets 
of  results.  The  shift  was  first  thought  to  be  possibly  due 
to  the  sparseness  of  the  mesh  used  in  the  models.  It  was 
however,  consistently  present  using  mesh  densities  of  up 
to  X/60. 

The  closeness  of  the  NEC  and  MININEC3  far- 
field  and  impedance  results  shown  here  suggests  that 
MININEC3  is  capable  of  the  comparable  emulation  of 
the  current  distribution  on  complex  structures  when 
coaq>ared  to  NEC.  The  codes  use  dissimilar  expansion 


and  testing  functions  to  model  the  current  distribution 
and  minor  differences  in  results  would  therefore  be 
expected.  The  level  of  agreement  obtained  though, 
especially  with  the  far-fields,  proves  the  validity  of 
using  MININEC3  to  model  complex  structures 
represented  as  wire  grid  models. 

6  Comparison  with  Experimental  Data 

As  a  final  test  of  the  ability  of  MININEC3  to 
model  complex  antenna  structures  accurately  we  used  an 
appropriate  set  of  experimental  data  as  a  reference. 
Bhattacharya  et  al  (1987)  measured  the  input  inqiedance 
of  a  monopole  antenna  mounted  on  a  conducting  metal 
box.  The  structure  was  a  five  sided  10cm  cube  attached 
centrally  to  a  lOScm  square  conducting  ground  plane. 
A  base-fed  6cm  vertical  monopole  was  mounted  at 
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various  podtioiis  oo  the  upper  fKe  of  the  cube.  The 
radiatioo  patUans  of  the  same  configuratioa  were  later 
measured  by  Chu  et  al  (1990)  in  the  same  frequency 
range  as  Bhattacharya,  l-2GHz. 

To  model  this  configuration  with  MININEC3 
a  wire  grid  modd  was  again  set  up.  Each  surface  of  the 
box  was  divided  into  a  6x6  meA  and  the  equal-area 
criterirm  was  enforced  to  calculate  the  wire  radius.  The 
resulting  wire  grid  model  is  shown  in  figure  14  with  the 
monopole  mounted  in  the  centre.  Two  simplificatiins 
w«e  also  made  with  the  MININEC3  model.  Firstly  the 
wires  of  the  grid  were  assumed  to  be  pofectly 
conducting  and  secondly  the  finite  conducting  ground 
plane  was  assumed  large  enou^,  in  term  of 
wavelengths,  to  be  replaced  by  a  infinite,  perfectly 
conducting  ground  plane. 


F%ure  14.  Wire  grid  representation  of  the 
conducting  box  with  centre-mounted 
vertical  monopole. 

Figures  IS  and  16  compare  the  MININEC3 
predicted  input  conductance  and  suscq>tance  with  the 
experimental  results  of  Bhattacharya  er  al.  It  is 
noticeable  that  both  the  computed  and  experimental 
results  have  a  similar  trmd  with  a  clear,  well  defined 
resonance.  For  both  G  and  B  however  the  peak  values 
differ  somewhat  and  there  is  a  noticeable  differmce 
betwemi  predicted  and  measured  resonant  frequency. 


FrtquMemy  (GHt) 

Figure  15.  Comparison  of  experimental  and 
MININEC3-predicted  input 
conductance. 


Figure  16.  Comparison  of  experimental  and 
MINlNEC3-predicted  input  susceptaiKe. 

To  test  if  this  difference  was  due  to  the 
inability  of  the  chosen  wire  grid  to  represent  the 
continuous  surface  accurately  two  further  models  were 
developed.  Firstly  the  mesh  density  was  doubled,  using 
a  12x12  grid  to  model  each  surface  of  the  cube. 
Seccmdly  the  twice  surface  area  rule  was  enforced  and 
this  yielded  a  larger  wire  grid  radius.  Both  of  these 
more  refined  models  however  only  marginally  improved 
the  agreement  of  the  MININEC3  results.  These  results 
therefore  proved  that  the  MININEC3  model  had 
converged  in  terms  of  the  necessary  grid  doisity  and 
area  foctor.  The  discrqrancies  are  therefore  probably 
due  to  two  further  factors.  Firstly,  the  assumption  made 
regarding  the  replacement  of  the  finite  ground  plane 
with  an  infinite,  perfectly  conducting  (me  may  not  be 
valid.  Secondly  external,  unspecified  factors  due  to  the 
measuremoit  system  are  another  probable  cause.  Cox 
(1991)  compared  the  measured  and  NEC-predicted  input 
impedances  of  a  number  of  HF  antennas  on  aircraft.  He 
attributed  differmces  which  existed  between  predicted 
and  experimental  input  inqredance  to  dissimilarities  in 


64 


the  modelled  and  actual  feed-points.  He  also  showed 
how  slight  differences  in  the  position  of  the  acttial  point 
of  measurement  introduced  a  predominantly  reactive 
differential  between  results.  Considering  the 
MININEC3  model  of  the  centre-mounted  monopole  the 
difference  between  the  predicted  and  experimental 
reactance  at  each  frequency  corresponded,  within  ±5%, 
to  an  inductive  offset  of  L  =  2.938nH.  Introducing  this 
inductance  in  series  with  the  MININEC3  predicted 
results  produced  the  G  and  B  plots  shown  in  figures  17 
and  18.  Clearly  this  empirically  determined  reactive 
shift  reduces  the  frequency  offset  of  the  resonance 
significantly.  Also  the  difference  in  the  magnitude  of  G 
is  reduced  substantially. 
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Figure  17.  Comparison  of  experimental  and 
MININEC3-predicted  input  conductance 
with  an  empirically  determined  inductive 
offset  in  the  predicted  results. 
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Figure  18.  Comparison  of  experimental  and 
MININEC3-predicted  input  susceptance 
with  an  empirically  determined  inductive 
offset  in  the  predicted  results. 


7  Conclusions 

This  paper  has  addressed  the  validation  of  an 
expanded  MININEC3  method  of  moments  computer 
code.  The  antenna  structures  used  in  the  validation  were 
particulary  appropriate  for  assessing  the  code's  ability 
to  model  complex  radiating  systems  accurately.  Until  its 
recent  conversion  from  BASIC  to  FORTRAN  it  was 
impossible  to  use  this  code  for  such  applications 
because  of  the  computer  CPU  time  and  compiler 
constraints.  The  three  techniques  used  in  this  validation 
exercise  show  that  MIN1NEC3  may  be  used  with 
confidence  to  model  electrically  large  structures. 
Firstly,  excellent  agreement  was  obtained  with 
analytical  results  for  the  input  impedance  of  a  V-dipole 
antenna.  This  was  achieved  even  at  the  most  acute  angle 
considered  of  30°,  provided  adequate  segmentation  was 
employed  at  the  junction  of  the  wires.  Secondly,  good 
agreement  was  obtained  with  the  NEC-predicted  results 
for  the  radiation  patterns  and  input  impedances  of  two 
wire  grid  models  of  typical  vehicular  type  structures. 
An  important  finding  however  was  the  noticeable 
frequency  shift  between  the  two  codes.  Although  this 
was  most  noticeable  because  of  the  near-resonant 
structures  that  were  used  in  the  simulation  it  must 
clearly  be  considered  when  employing  the  code.  Finally 
the  good  agreement  with  the  experimental  data  for  the 
input  impedance  of  a  monopoie  antenna  on  a  conducting 
box  again  demonstrated  the  usefulness  of  this  modified 
MININEC3  code. 
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ABSTRACT 

Examples  are  presented  of  the  application  of  a  moment-method  modelling  procedure  (using 
NEC)  to  calculate  the  near-field  dependent  parameters  of  aircraft  HF  antennas.  The  overall 
aim  of  the  work  from  which  the  examples  are  drawn  is  to  provide  a  ‘recipe’;  the  application 
of  which  will  generate  predictions  of  parameters  of  interest  which  are  both  credible  and 
sufficiently  accurate  for  engineering  purposes.  A  feature  of  the  procedure  is  that  no 
empirically  derived  information  is  incorporated.  The  examples  encompass  a  range  of 
generic  aircraft  and  antenna  types:  both  fixed  and  rotary  wing  airframes,  and  electric  (wire) 
and  magnetic  (loop  and  notch)  primary  radiators.  Particular  emphasis  is  placed  upon 
validation  of  the  predictions.  This  is  performed  by  direct  measurement  where  possible,  but 
where  this  is  not  practicable  simple,  but  independent,  calculations  are  performed  based  upon 
equivalent  circuits.  As  for  the  moment-method  procedure,  these  corroborative  calculations 
do  not  depend  upon  empirical  information. 

INTRODUCTION 

At  HF  (2-30  MHz)  airframes  are  resonant  at  numerous  frequencies  throughout  the  band  and 
modelling  procedures  must  be  sufficiently  sound  to  predict  such  resonances  and  their  effect 
upon  parameters  of  interest  Thus,  near-fields  must  be  calculated  accurately  for  this  general 
reason  as  well  as  the  more  particular  one  of  predicting  the  actual  value  of  parameters  which 
depend  upon  the  near-field  {eg  antenna  reactance,  coupling  between  antennas  etc).  Accurate 
calculation  of  the  near-field  requires  good  solutions  for  the  body  surface  currents  and  this 
implies  high  sampling  densities.  There  are,  however,  limitations  to  the  sampling  density 
which  can  be  contemplated.  The  more  obvious  practical  limitations  are: 

(i)  time  taken  to  construct  mathematical  models; 

(ii)  run  times; 

(iii)  factors  depending  upon  machine  precision,  basis  functions,  matrix  size,  and 
the  possible  consequences  for  solution  stability. 

For  those  involved  in  practical  calculations  pertinent  questions  therefore  are: 

(a)  How  high  does  the  sampling  density  need  to  be  for  engineering  purposes? 

(b)  How  cT^ble  are  the  results? 

The  answers  to  these  questions  depend  of  course  upon  particular  circumstances  and 
requirements,  but  it  is  hoped  that  the  following  examples  might  be  of  some  use  in  this 
regard. 

In  evolving  the  methods  described  here,  it  was  considered  important  that  the  modelling 
procedure  would  be  truly  predictive:  that  is,  it  should  not  involve  any  empirically  derived 
information.  For  the  examples  cited  below,  construction  of  the  mathematical  model  is 
performed  only  with  the  aid  of  general  arrangement  drawings  and  (when  available) 
photographs  or  models.  The  reason  for  this  is  that  predictions  of  the  likely  performance  of  a 
proposed  antenna  system  are  usually  required  early  in  a  project  cycle  prior  to  the  airframe 
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being  available  for  measurement,  or  even  in  existence  at  all.  Indeed  if  circumstances  are 
otherwise  it  is  sometimes  better  to  generate  the  required  information  by  other  means.  Thus, 
the  same  modelling  procedure  was  used  without  ad  hoc  modification  for  each  of  the 
examples,  and  is  as  follows. 

MOMENT-METHOD  MODELLING  PROCEDURE 

For  aircraft-like  (that  is,  non-voluminous)  structures  in  the  HF  band,  the  solution  stability 
problems  of  the  magnetic  field  integral  equation  (MFIE)  preclude  its  use  unless  prohibitively 
high  sampling  densities  are  adopted.  Use  of  the  electric  field  integral  equation  was  therefore 
consider^  necessary.  The  aircraft  surface  is  modelled  as  a  wire-mesh  and  discretized  with 
the  aim  of  producing  an  approximately  square  mesh  (becoming  triangular  when  required  by 
Itxal  topology)  with  a  mesh  side  of  atout  0.0025  wavelength  (0.3  to  0.4  m)  at  the  bottom  of 
the  band  (2  MHz).  The  solution  code  is  a  single  precision  version  of  NEC3(1)  with 
modifications  to  facilitate  running  on  a  Cray  1  machine.  The  sampling  density  was  made 
one  match  point  per  mesh  side.  The  radii  of  the  mesh  wires  are  such  that  the  resulting 
surface  area  of  each  wire  when  notionally  ‘unwrapped’  is  the  same  as  that  of  the  surface  that 
it  replaces.  In  effect  this  represents  an  attempt  to  satisfy  (approximately)  the  MFEE,  the  first 
term  of  which  represents  Ampere’s  law  at  a  portion  of  a  perfectly  conducting  closed  surface. 
That  is,  the  magnetic  field  boundary  conditions  are  approximately  satisfied,  in  addition  to 
the  electric  field  boundary  conditions  that  are  imposed  at  the  match  points. 

VALIDATION  OF  MOMENT-METHOD  PREDICTIONS 

A  means  of  quantitatively  assessing  the  vahdity  of  moment-method  calculations  for  small- 
scale  problems  is  to  examine  the  convergence  of  solution.  That  is,  the  discretization  is 
successively  refined  and  the  predicted  result  of  some  parameter  of  interest  examined  after 
each  calculation.  The  anticipation  is  that  the  result  at  each  stage  will  converge  asymptotically 
to  the  ‘true’  value.  This  is  not  a  practical  proposition  for  large  scale  problems  of  the  type 
considered  here:  the  constraints  of  successively  longer  periods  for  model  construction  and 
run  times  preclude  this.  Thus,  experimental  and/or  independent  theoretical  means  must  be 
sought  for  validation. 

A  major  problem  in  assessing  the  credibility  of  predictions  concerning  aircraft  HF  antennas 
by  experimental  means  is  the  frequently  encountered  impracticability  of  measuring  some 
parameters  of  engineering  interest.  An  obvious  example  is  radiation  resistance:  at  the  low 
frequency  end  of  the  band  this  is,  in  the  vast  majority  of  cases,  much  smaller  than  the 
airframe  loss  resistance  and  cannot  be  determined  by  a  measurement  made  at  the  antenna 
terminals.  It  is  clearly  not  a  practical  proposition  to  measure  the  total  radiated  power  in 
order  to  deduce  the  radiation  resistance.  In  addition,  to  measure  those  parameters  which 
could  be  determined  experimentally  is  often  prohibitively  expensive  because  of  engineering 
woric  which  must  first  be  performed  upon  the  aircraft.  This  is  particularly  true  of  smaller 
aircraft  where  antenna  terminals  are  often  located  at  positions  which  are  inaccessible  when 
the  aircraft  is  in  flight.  Thus,  independent  supporting  evidence  against  which  the  credibility 
of  predictions  can  be  assessed  is  usually  meagre  and  so,  whatever  direct  or  circumstantid 
information  and  physical  insight  which  is  available  must  be  exploited:  these  difficulties  have 
been  considered  elsewhere  (2). 

In  the  examples  considered  below  only  in  one  case  was  direct  (that  is,  in  flight)  experimental 
evidence  obtainable  for  corroboration  of  near-field  parameters.  This  was  because  the 
subject  was  a  medium-size  passenger  aircraft  and  the  terminals  were  easily  accessible  from 
within  the  fuselage.  In  another  example,  terminal  measurements  were  only  practicable  for 
the  aircraft  locat^  on  the  ground;  this  involves  additional  effects  and  complications  which 
are  not  easy  to  quantify.  In  another  case  indirect  evidence  is  sought  from  radiation  patterns: 
predicted  patterns  were  compared  with  those  derived  from  measurements  made  upon  a 
physical  scale  model,  and  upon  a  real  aircraft  in  flight.  Although  a  far-field  dependent 
property,  radiation  patterns  show  evidence  of  airframe  resonances  and  therefore  give  an 
indication  of  the  accuracy  with  which  near-fields  are  calculated. 
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Where  experimental  corroboration  is  not  practicable,  ^uivalent  circuit  models  are  derived  to 
assess  the  credibility  of  the  moment-method  calculations.  The  circuit  parameters  for  these 
models  are  derived  in  a  logical  manner  and  not  simply  adjusted  to  fit  the  moment  method 
predictions.  These  circuits,  however,  involve  the  application  of  insight  in  their 
construction,  considerable  simplification  to  maintain  tractability  (and  therefore  the 
introduction  of  an  element  of  subjectivity),  and  can  only  be  convincingly  produced  for 
certain  antenna  types  and  configurations.  It  is  of  course  unreasonable  to  expect  good 
agreement  between  the  moment-method  calculations  and  those  derived  from  an  equivalent 
circuit.  In  fact  the  circuit  nxxiel  can  really  serve  to  give  only  an  indication  as  to  whether  the 
moment-method  calculations  are,  or  are  not,  credible.  One  reason  for  this  is  that  in  those 
cases  where  the  construction  of  a  simple  equivalent  circuit  is  possible,  the  antennas  are 
located  near,  but  not  precisely  at,  the  electrical  centre  of  the  fuselage.  The  readily  obtainable 
usable  expressions  or  tabulated  values  for  Z ,  the  impedance  of  ^e  dipole  representing  the 
fuselage  (see  below)  are,  however,  derived  for  centre  excitation.  Fortuitously  for  the  two 
examples  (2,9)  in  which  equivalent  circuits  are  constructed  the  antennas  are  located  at 
approximately  40%  along  the  electrical  length  of  the  fuselage.  In  addition,  in  both  of  these 
cases  the  resulting  circuit  component  values  are  such  that  Z  does  not  exercise  a  dominating 
influence  on  the  final  result.  Thus,  notwithstanding  these  shortcomings,  it  is  expected  that, 
when  applicable,  quantitative  agreement  between  the  moment-method  calculation  and  the 
equivalent  circuit  would  be  to  better  than  an  order  of  magnitude.  By  this  means  therefore,  a 
relatively  crude,  but  independent,  piece  of  theoretical  supporting  evidence  can  be  produced. 

The  calculations  (that  is,  moment-method  and  equivalent  circuit)  share  a  common  deficiency: 
neither  incorporate  any  allowance  for  the  effects  of  non-zero  impedance  of  the  mechanical 
joints  which  unavoidably  exist  on  real  aircraft.  As  far  as  HF  airci^t  antenna  performance  in 
general  is  concerned  these  extraneous  impedances  are  a  perennial  source  of  difficulty.  This 
is  because  throughout  much  of  the  band  most  aircraft  are  not  electrically  large.  This 
inevitably  implies  that  radiation  resistance  and,  in  practice  radiation  efficiency,  tends  to  be 
relatively  small:  estimated  values  of  below  0.1%  are  not  uncommon  at  the  bottom  of  the 
band.  Not  only  are  joint  impedances  unknown  (and,  for  practical  purposes,  unknowable) 
but  they  vary  between  nominally  identical  aircraft  and,  in  the  mechanically  harsh 
environment  experienced  by  aircraft,  change  over  time.  Thus  the  assumption  of  perfect 
conductivity  for  calculations  is  practically  inevitable  as  far  as  such  joints  are  concerned. 
Since  the  extraneous  resistance  due  to  joints  is  usually  dominant  over  much  of  the  band 
(experience  indicates  that  measured  terminal  resistance  is  much  greater  than  that  calculated 
on  the  assumption  of  no  joints  at  all),  there  is  little  to  be  gained  by  assuming  other  than  a 
perfectly  conducting  airframe.  That  is,  the  calculated  terminal  resistance  is  the  radiation 
resistance  here. 

In  one  of  the  airframe  examples  there  are  significant  portions  of  the  aircraft  skin  made  from 
non-metallic  material.  It  was  noted,  however,  that  large  areas  of  metal  existed  immediately 
behind  the  skin.  It  was  therefore  considered  realistic,  at  least  at  HF,  to  neglect  the  skin  and 
treat  the  surface  as  perfectly  conducting  and  electrically  bonded  to  the  rest  of  the  airframe.  It 
is  expected  that  such  bonding  of  internal  metal  structures  will  be  arranged  to  afford 
protection  against  lightning  strike. 

Example  1:  Terminal  reactance  of  a  medium-size  passenger  aircraft. 

Three  ‘long-wire’  antennas  are  fitted  to  the  aircraft:  the  arrangement  is  shown  in  Fig  1.  The 
model  comprises  7736  segments  and  this  is  the  number  of  uidmown  currents.  This  antenna 
type  and  configuration  presents  potential  sources  of  computational  difficulty.  These  are 
considered  more  fully  elsewhere  (3)  but  briefly  they  are  connected  with: 

( 1 )  The  very  close  electrical  proximity  of  the  match  points  in  the  region  of  the  antenna 
feed  point  and  near  the  antenna  anchoring  point  at  the  tail  fin.  The  match  point  separation 
distances  are  around  1/3000  and  1/1500  of  a  wavelength  respectively  for  these  two  regions. 


69 


(2)  The  problem  of  distribution  of  electric  charge  between  wires  of  disparate  radii  at 
junctions.  An  incidental  advantage  of  adopting  a  very  fine  mesh  is  that  the  radii  of  the  wires 
simulating  the  fuselage  become  smaller.  In  addition,  the  ratio  of  the  radius  of  the  wires 
simulating  the  fuselage  to  that  of  the  antenna  wire  falls  as  the  mesh  is  made  finer.  For  the 
junction  treatment  (based  upon  a  study  of  the  tapered  antenna  by  Wu  and  King  (Ref  4))  used 
in  NEC,  when  the  wires  become  electrically  very  thin  the  linear  charge  density  at  each  of  the 
wires  at  a  junction  is  the  same  and  does  not  depend  strongly  on  the  wire  radii.  Thus  for 
finer  meshes  it  is  anticipated  that  the  solutions  become  less  dependent  upon  the  arrangement 
made  for  dealing  with  junctions.  Even  with  the  high  sampling  densities  used  here, 
however,  the  ratio  referred  to  above  is  still  about  19.  The  concern  is  that  if  unrealistic 
charge  distributions  at  the  base  (that  is,  the  terminals)  of  a  wire  antenna  are  computed,  these 
may  be  expected  to  have  a  significant  effect  upon  the  calculated  terminal  impedance. 

Comparison  between  prediction  and  measurement  for  the  terminal  reactance  is  depicted  in 
Fig  2  with  allowance  made  for  the  presence  of  the  lightning  spark  gap  and  other  measured 
stray  reactances  present  in  the  real  installation  (in  effect  a  ‘front-end’  network  to  the  actual 
antenna  terminals)  but  which  cannot  easily  be  incorporated  directly  in  a  moment-method 
model  (see  Ref  3).  It  is  seen  that  the  resonances  are  predicted  to  within  about  10%  of  the 
measured  values. 

Example  2:  Calculation  of  maximum  mutual  coupling  between  two  HF  loop  antennas 
mounted  upon  a  helicopter  at  the  bottom  of  the  HF  band. 

The  arrangement  is  depicted  in  Fig  3.  There  are  2250  segments  in  the  model.  The 
maximum  mutual  coupling  is  a  function  of  the  Y  parameters  of  the  configuration  (see  Ref  2) 
and  its  estimation  is  sensitive  to  the  accuracy  with  which  these  parameters  are  calculated. 
Validation  of  the  moment-method  calculations  was  attempted  by  the  derivation  of  an 
equivalent  circuit.  Below  the  first  airframe  resonance  physical  considerations  suggest  that 
the  airframe  resembles  an  electric  dipole  which  is  shunt  fed  by  the  driven  loop  antenna. 
Thus  the  dominating  elementary  radiating  modes  would  be  expected  to  be  as  indicated  in 
Fig  4.  In  this  spectral  region,  and  for  this  type  of  antenna,  it  is  expected  that  the  coupling 
mechanism  between  the  driven  antenna  and  the  non-driven  antenna  is  primarily  magnetic. 
That  is  the  driven  antenna  excites  the  other  principally  by  transformer  action.  The  equivalent 
circuit  for  the  arrangement  is  illustrated  in  Fig  5.  An  approximate  estimate  of  the  mutual 
inductance,  M  ,  between  the  two  tail-cone  mounted  antennas,  was  made  by  treating  the 
arrangement  as  a  DC  quasi-two-dimensional  problem  and  performing  a  moment-method 
calculation.  An  elementary  calculation  of  the  remaining  inductances  was  made  using  the 
method  of  images  in  conjunction  with  a  procedure  which  has  been  termed  by  Roters  (5)  as 
‘estimating  the  permeances  of  probable  flux  paths’.  The  radiation  resistance,  rj,  of  the  loop 
is  estimated  using  the  well-known  expression  for  an  electrically  small  magnetic  dipole  with 
allowance  made  for  the  fact  that  it  is  mounted  on  a  conducting  surface.  A  more  detailed 
account  of  the  derivation  of  the  equivalent  circuit  is  given  in  Ref  2.  The  impedance,  Z ,  of 
the  dipole  representing  the  fuselage  was  estimated  using  the  following  analytical  expression 
(6). 


Z  =  R(k«)-j[120(ln(£/a)-l)cot(k£)-X(k£)],  (1) 

where  £  is  the  half  length  of  the  dipole,  a  its  radius,  and  k  the  free  space  wave  number 
of  the  excitation,  the  functions  R(k£)and  X(k^)  are  tabulated  in  Ref  (6).  This  expression 
is  not,  however,  valid  at  frequencies  higher  than  the  first  resonance. 

Comparison  of  the  maximum  mutual  coupling  calculated  by  the  two  methods  is  shown  in 
Fig  6.  The  error  bars  in  Fig  6  are  a  consequence  of  the  uncertainty  in  assigning  a  figure  to 
the  equivalent  radius  of  the  dipole  representing  the  fuselage,  and  indicate  the  extreme 
possible  values. 

An  indicator  of  solution  quality  can  also  be  made  by  examining  compliance  with  reciprocity. 
This  is  shown  in  the  following  table. 
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Frequency 

(MHz) 

Current  at  antenna  2  terminals 
while  short-circuited  and 
antenna  1  driven  with  excitation 

Current  at  antenna  1  terminals 
while  shOTt-circuited  with 
antenna  2  driven  with  excitation 

of  IV  (A) 

of  IV  (A) 

2 

-0.72X  10-7+ j  0.15  X  10-3 

-0.85X  10-7+ j  0.17  X  10-3 

3 

0.59X  10-6  + j  0.12  X  10-3 

0.62X  10-6+ j  0.13  X  10-3 

4 

0.57X  10-5+ j  0.11  X  10-3 

0.61  X  10-5 +j  0.12  X  10-3 

5 

0.42  X  10-4+ j  0.13  X  10-3 

0.46X  10^+j0.14x  10-3 

6 

0.51  X  10-4 +j  0.33  X  10-4 

0.73X  10^+j0.42x  10^ 

Table  1  Compliance  with  reciprocity 

Example  3:  Determination  of  antenna  impedance  for  a  notch  antenna  installed  on  a 
variable  geometry  aircraft. 

The  model  is  constructed  in  such  a  way  that  it  can  be  reconfigured  as  a  real  aircraft  of  this 
type  might  change  its  geometry.  That  is,  alterations  to  the  model  for  various  positions  of  the 
moveable  parts  correspond  exactly  to  the  geometrical  changes  of  the  wing  or  taileron 
positions  alone;  there  is  no  remeshing  or  any  other  ad  hoc  alternations  made  to  the  model 
(9,10).  Figs  7  and  8  indicate  the  model  with  the  wings  in  the  fully  forward  and  swept  back 
positions.  There  are  3293  segments  in  the  model.  It  was  considered  that  the  following 
factors  might  be  a  potential  source  of  computational  difficulty. 

(i)  The  moveable  surfaces  clear  the  rest  of  the  airframe  by  electrically  ve^  small 
distances  (of  the  order  of  1/1500  of  a  wavelength  at  2  MHz).  In  addition  this 
necessitates  the  use  of  very  short  segments  (representing  pins  for  mechanical 
pivoting)  in  the  model  which  are  of  a  similar  electrical  length  to  this  figure. 

(ii)  The  area  of  the  notch  antenna  is  of  the  order  of  the  mesh  size.  In  the  context 
of  running  NEC  on  32-bit  machines,  difficulties  have  been  indicated  for  problems 
involving  small  loops  (7).  Although  the  work  reported  here  was  performed  with 
single  precision  on  64-bit  machines,  it  is  considered  that  calculations  should  be 
regarded  with  caution  particularly  in  view  of  the  large  and  complex  nature  of  the 
structures. 

Validation  of  predictions  of  the  terminal  impedance  made  for  the  aircraft  in  flight  was 
attempted  by  two  methods.  The  first  was  to  derive  an  equivalent  circuit  in  a  manner  similar 
to  that  indicated  for  the  previous  example.  The  second  was  by  measurement  of  terminal 
impedance  with  the  aircrtrft  on  the  ground. 

The  equivalent  circuit  is  show  in  Fig  9.  As  indicated  above,  the  use  of  equation  (1)  to 
calculate  Z  is  limited  to  frequencies  below  the  first  airframe  resonance.  In  addition  to  using 
equation  (1),  Z  was  also  estimated  assuming  that  the  effect  of  dipole  thickness  on  the 
antenna  terminal  impedance  is  not  too  great.  This  is,  of  course,  a  rather  gross 
approximation,  but  it  does  permit  values  of  Z  to  be  obtained  beyond  the  first  airframe 
resonance  using  tabulated  results  using  Wu’s  theory  (8). 

The  ground  measurements  were  performed  using  a  vector  impedance  meter  (VIM).  Such 
measurements  are  not  entirely  straightforward  because  the  effects  of  the  cable  which  is 
required  to  connect  the  VIM  to  the  antenna  terminals  must  be  removed  from  the  data  in  order 
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to  extract  the  terminal  impedance.  This  aspect  is  reported  in  more  detail,  along  with  checks 
of  the  experimental  method,  in  Ref  10. 

Comparison  between  the  NEC  calculation,  the  equivalent  circuit  model,  and  the  ground 
measurements  are  depiaed  in  Fig  10  for  the  terminal  reactance.  As  indicated  above  it  is  not 
possible  to  extract  the  radiation  resistance  from  the  measured  terminal  resistance;  hence  for 
the  former  parameter,  the  moment-method  predictions  can  only  be  compared  with  those  of 
the  equivalent  circuit  (see  Fig  11).  Elementary  considerations  suggest  that  the  reactance 
should  increase  monotonically  throughout  the  band  and  remain  positive.  This  is  observed 
for  the  calculations  and  measurements.  It  is  seen  that  the  values  derived  by  measurement  are 
about  36%  below  those  of  the  NEC  calculation  over  that  part  of  the  band  which  is 
approximately  linear.  This  figure  is  comparable  with  a  similar  difference  (26%)  between 
measured  and  theoretical  values  (from  an  analytical  expression)  performed  as  a  check  of  the 
measurement  procedure  and  made  on  a  circular  loop  having  an  area  comparable  to  that  of  the 
notch  (10). 

It  is  expected  from  simple  physical  considerations  that  the  longitudinal  airframe  resonances 
will  occur  around  those  frequencies  at  which  the  electrical  length  of  the  airframe  is  an  integer 
multiple  of  half  the  radiation  wavelength.  They  are  approximately  7.6, 15.2  and  22.8  MHz 
for  this  example.  The  form  of  the  equivalent  circuit  (see  Fig  9)  and  the  fact  that  the  radiation 
resistance  of  the  dipole  is  always  significantly  greater  than  that  of  the  notch  (ri),  suggests 
that  such  effects  are  expected  to  be  noticeable  in  the  NEC  calculations  even  though 
resonance  is  a  phenomenon  associated  with  the  near-field  whereas  radiation  resistance  is  a 
far-field  property.  It  is  seen  in  Fig  1 1  that  the  results  of  the  NEC  calculation  show 
disturbances  which  occur  close  to  the  above  values  of  resonant  frequency.  Such  effects  re 
not  observed  for  the  reactance  (Fig  10)  because  the  component  values  are  such  that  the 
shunting  effect  of  Z  is  not  such  as  to  be  large  at  any  frequency  in  the  band;  not  even  at  a 
series  resonance.  The  pronounced  oscillations  seen  in  Fig  1 1  for  the  equivalent  circuit 
model  are  a  result  of  the  thin  wire  approximation  for  the  fuselage.  Such  oscillations  would 
not  be  expected  to  be  so  marked  in  the  case  of  a  large  diameter  dipole  such  as  that  forming 
the  real  fuselage  and  it  is  seen  in  Fig  1 1  that  they  are  only  just  perceptible  on  the  curve  for 
the  NEC  calculation.  Interestingly,  the  two  methods  display  a  reasonable  agreement  at  the 
series  resonances  (7.6  and  22  NMz).  At  such  frequencies  the  representation  of  the  fuselage 
by  a  thin  wire  would  not  be  expected  to  exert  a  great  influence  to  the  radiation  resistance 
estimated  using  the  equivalent  circuit. 

Example  4;  Inference  of  first  longitudinal  airframe  resonant  frequency  from  radiation 
patterns. 

This  helicopter  was  fitted  with  a  loop  antenna  and  is  as  depicted  in  Fig  12.  There  are  2682 
segments  in  the  model.  The  example  is  included  to  indicate  that,  as  far  as  validation  is 
concerned,  radiation  patterns  can  be  used  to  indicate  quality  of  solution  even  though 
radiation  patterns  depend  upon  the  far-field.  This  is  because,  at  HF,  airframes  are  generally 
resonant  at  numerous  frequencies  throughout  the  band.  For  this  example  measurements  of 
radiation  patterns  were  made  both  using  a  scale  physical  model  and  with  a  real  aircraft  in 
flight.  Details  of  this  study  along  with  a  set  of  radiation  patterns  throughout  the  HF  band  is 
recorded  elsewhere  (11).  The  first  airframe  resonance  is  usually  apparent  in  azimuthal 
radiation  patterns  for  horizontal  polarisation.  It  is  characterised  by  a  deep  null  in  the  general 
nose  and  tail  directions.  The  nulls  do  not  lie  exactly  on  the  nose-tail  axis  unless  the  antenna 
itself  is  also  located  symmetrically  about  this  line.  Greater  complexity  of  the  patterns  at 
higher  frequency  resonances  renders  such  unambiguous  interpretation  less  likely  than  in  the 
case  of  the  first  airframe  resonance.  This  is  due  to  the  corresponding  increase  in  the  number 
of  possible  radiation  modes  which  will,  in  general,  be  present.  For  the  antenna  orientation 
adopted  here,  the  length  of  the  probable  RF  current  path  indicates  that  the  first  longitudinal 
airframe  resonance  should  be  in  the  region  of  6  MHz.  Figs  13-15  indicate  a  succession  of 
such  patterns  as  the  frequency  is  increased  through  the  first  airframe  resonance.  These 
figures  show  the  anticipated  resonance  occurs  around  5.6  MHz.  In  passing  it  is  noted  that 
the  patterns  display  the  fact  that  the  antenna  is  located  on  the  left  side  of  the  aircraft.  This  is 
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seen  from  the  slightly  higher  levels  of  radiation  on  this  side  and  if  no  experimental  evidence 
was  available  would  provide  a  simple  but  nevertheless  useful  additional  indicator  of  the 
credibility  of  the  calculations. 

Examination  of  Figs  13-15  indicate  in  addition  an  aspect  which  is  of  some  practical 
importance;  namely  that  there  is  a  level  of  precision  in  calculations  beyond  which  there  is 
only  limited  practical  benefit  in  proceeding  at  least  in  the  case  of  radiation  pattern  prediction 
at  HF.  This  is  due  to  the  difficulty  of  m^ng  measurements  reliably  and  repeatably  with 
real  aircraft  and  is  in  pan  due  to  the  stochastic  nature  of  the  propagation  conditions  in  this 
spectral  region.  Thus,  there  is  some  uncertainty  as  to  what  the  ‘correct  answer’  actually  is. 

CONCLUSIONS 

For  the  examples  presented,  the  moment-method  calculations  seem  generally  consistent  with 
such  attempts  at  validation  as  could  be  made,  and  it  is  consider^  that  they  enhance  the 
degree  of  confidence  which  can  be  placed  in  the  prediction  of  near-field  dependent 
parameters  for  large  scale  problems  involving  aircraft.  For  the  aircraft  HF  case  at  least,  the 
difficulty  of  quantifying  the  electrical  properties  of  joints  present  on  the  airframe,  and  the 
variability  of  the  propagation  conditions,  limit  the  precision  which  is,  in  practice,  useful 
when  calculations  are  required  for  the  estimation  of  certain  system  parameters  such  as  link 
budgets  and  for  the  prediction  of  radiation  patterns. 
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Abstract 

A  direct  method  for  the  numerical  solution  of 
Maxwell's  equations  in  the  frequency  domain 
(sinusoidal  steady  state),  for  radiation  problems  of 
cylindrical  symmetry,  is  incorporated  into  a  finite 
difference  program  called  HERZ.  coded  in  FORTRAN. 
The  termination  of  the  fields  on  the  outer  boundary  of 
the  problem  domain  is  accomplished  via  an  absorptive 
layer,  labelled  a  stealth  layer,  which  attenuates  the 
incident  fields  to  insignificant  amounts  in  a  small 
computing  space,  without  causing  significant 
reflections  that  would  disturb  the  near  field  solutions. 
No  radiation  boundary  condition  is  specified  on  the 
problem  domain,  rather,  the  transmission  of  the  fields 
to  the  stealth  layer  is  optimized  by  specifying  its 
electrical  properties.  To  validate  HERZ  as  a  useful 
electromagnetic  modeller  for  thin  wire  antennas, 
driving  point  admittances  and  surface  current 
distributions  for  several  configurations  were  compared 
to  both  theoretical  and  measured  values,  with  good 
agreement. 

Introduction 

In  the  numerical  modeling  of  radiation  problems, 
particularly  those  of  isolated  thin  wire  antennas, 
specific  radiation  boundary  conditions  allow  finite 
difference  or  finite  element  techniques  to  be  applied 
(Ramahieta/.  1991).  The  implementation  of  non-local 
boundary  conditions  such  as  boundary  integral 
formulations  or  moment  method  techniques  is  useful 
for  infinite  homogeneous  configurations  but  bounded 
methods  are  typically  incorporated  in  conjunction  with 
these  techniques  if  inhomogeneities  exist,  such  as  in 
the  case  of  a  dielectrically  coated  antenna  (McDonald 
and  Wexler  1972,  Morgan  et  al.  1984,  Yuan  et  al. 
1990).  To  avoid  the  complication  of  coupling  a 
bounded  problem  to  an  unbounded  problem,  specific 
radiation  boundary  conditions  are  formulated  in 
conjunction  with  finite  element  techniques  (Sumbar  et 
al.  1990). 

While  finite  element  techniques  allow  for  very  flexible 
grid  density  and  are  straightforward  in  their  treatment 
of  inhomogeneites,  the  complex  nodal  grids  which 
must  be  generated  lack  the  physical  transparency  and 


simplicity  of  a  finite  difference  grid  system.  In  this 
paper  an  absortiing  layer  adjacent  to  the  problem  space 
boundary  is  used  in  place  of  an  explicit  radiation 
boundary  condition.  Ratho’  than  specifying  a  radiation 
condition  to  t^minate  the  fields  of  the  antenna,  shown 
in  Fig.  1,  on  the  problem  domain,  the  transmission  of 
the  fields  from  the  problem  domain  to  a  lossy  layer 
enveloping  the  problem  domain  is  modelled,  as 
illustrated  in  Fig.  2.  The  problem  domain's  boundaries 
are  then  in  an  interior  grid  location,  as  the  grid 
encompasses  both  the  problem  domain  and  it's 
surrounding  lossy  layer,  which  will  now  be  referred  to 
as  the  stealth  layer.  The  near-field  problem  can  then 
be  solved  at  the  expense  of  the  far-field  solutions, 
which  will  be  attenuated  in  the  stealth  layer,  without 
disturbing  the  field  distributions  within  the  initial 
problem  domain.  A  related  technique  for  truncating 
the  solution  domain  in  finite  element  scattering 
problems  has  been  discussed  in  the  literature  (Jin  et  al. 
1991, 1992). 

This  methodology  is  incorporated  into  a  FORTRAN 
program,  called  HERZ  (H  field  E  field,  R  -  Z 
geometry),  which  is  executed  on  a  Macintosh  II 
personal  computer.  The  finite  difference 
implementation  of  Maxwell's  equations  and  the  use  of 
the  absorptive  stealth  layer  are  described  in  the 
following  sections.  Several  antenna  configurations 
are  modelled  with  HERZ  to  validate  this  approach, 
comparing  field  related  values  such  as  driving  point 
admittances  and  surface  current  distributions  to 
published  theoretical  and  measured  results. 

Finite  difference  implementation  of 
Maxwell's  equations 

Cylindrically  symmetric  metallic  antennas  operating  in 
cylindrically  symmeoic  inhomogeneous  materials,  with 
geometries  as  shown  in  Fig.  1,  will  have  simplified 
electric  and  magnetic  field  orientations  when  the 
electrical  properties  of  these  materials  are  isotropic. 
The  electric  fields  will  have  radial  and  axial 
components  (r,z),  and  the  magnetic  fields  only  an 
azimuthal  component  (0). 
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The  general  time-harmonic  solutions  will  be  of  the 
form 

E  =  Er<r^)ar  +  Ez(r.z)aj  (1) 

H  =  H0(r,z)a0.  (2) 

Assuming  an  el^^  time  dependence,  the  integral  form 
of  Maxwell's  equations  can  be  expressed  as 

^Edl  =  -yw/x/jHdS  (3) 

and  ^H  dI  =  (a-i-y<»e)J|E  dS  .  (4) 

The  solution  technique  adopted  uses  these  integral 
equations  to  form  surfaces  and  contours  accommodated 
by  the  finite  difference  grid  structure  (Albani  and 
Bemardi  1974),  providing  an  algorithm  to  compute  the 
azimuthal  magnetic  field  component,  H0,  at  the  center 
of  every  gridblock.  (3)  was  discretized  to  express  the 
H0  field  in  a  gridblock  center  in  terms  of  the  four  E 
field  values  which  surround  it,  as  is  displayed  in  Fig.  3, 
the  finite  difference  grid  implementation  in  HERZ. 
The  resulting  linear  equation  is  of  the  form 

aij  Ezij  +  bij  Erij+i  -i-  Cij  Ezj+i  j  + 

dij  EriJ  =  eij  H0y  (5) 

where  the  coefficients  aij  to  eij  depend  on  the  grid 
dimensions  Ari  and  Azj  (Nachai  el  al.  1992).  The 
gridblocks  have  been  assumed  to  be  small  enough  that 
the  electric  field  components  are  constant  along  each 
line  segment  of  the  gridblock  ij  and  the  magnetic  field 
is  uniform  over  its  cross  section.  In  this  equation 

rH0  =  H0' 

where  r  is  the  radial  distance  to  the  center  of  the 
gridblock  where  H0  is  being  calculated. 

By  applying  (4)  to  appropriate  surfaces  and  contours, 
the  various  electric  field  components  in  (S)  can  be 
written  in  terms  of  the  magnetic  fields  at  the  centers  of 
the  surrounding  gridblocks.  For  example,  with 
reference  to  Fig.  3,  E^i  j  can  be  expressed  in  terms  of 
H0'i.ij  and  H0'ij  by  applying  (4)  to  the  contour 
consisting  of  two  concentric  circular  paths  of  radii  ri.] 
and  n  to  compute  the  total  axial  current  (conduction 
and  displacement)  through  the  enclosed  surface.  The 
equation  that  results  for  the  magnetic  field  is 

AiJ  H0'ij  +  Bij  H0'i.i  j  +  Cij  H0'i+i  j  + 

I^i j  H0'ij-i  +  Eij  H0'ij+i  =  0  (6) 


where  the  coefficients  Aj  j  to  Eij  are  given  in  terms  of 
the  frequency  o),  the  electrical  properties  a,  e  and  and 
the  radial  locations  and  dimensions  of  the  pertinent 
gridblocks  (Nachai  et  al.  1992).  The  discretizations  of 
(3)  and  (4)  are  particularly  transparent  and  easy  to 
implement  because  all  electric  and  magnetic  field 
components  are  tangential  to  material  interfaces  and 
the  difficulties  that  can  be  encountered  in  handling 
inhomogeneous  problem  domains  in  three  dimensions 
do  not  arise. 

Each  of  the  N  gridblocks,  referred  to  by  indices  iJ,  will 
have  an  associated  equation  of  the  form  of  (6), 
providing  the  N  unknown  H  fields  with  N  linearly 
independent  equations,  resulting  in  a  determined 
system.  To  completely  specify  the  problem,  conditions 
must  be  provided  on  the  problem  boundary.  Excitation 
by  electric  fields  on  the  grid's  outer  domain,  the 
Neumann  boundary  condition,  requires  an  adjustment 
of  the  LHS  of  (6),  as  one  of  the  terms  Er  or  Ez  in  (5)  is 
then  predetermined.  The  excitation  electric  field 
appears  on  the  RHS  of  (6). 

Electric  field  excitation  is  specified  exclusively  on  the 
perimeter  of  the  problem  domain,  where  it  may,  for 
instance,  be  set  to  zero  to  signify  the  presence  of  a 
perfectly  conducting  boundary.  At  all  other  locations 
in  the  problem  domain  the  electric  field  excitation  is 
zero.  Magnetic  field  excitation  may  be  applied  at  the 
center  of  any  gridblock  by  specifying  a  value  of 

The  set  of  N  equations  for  H0'  represents  a  banded 
symmetric  system  which  can  be  stored  in  the  computer 
in  a  compressed  coefficient  mauix,  K.  The  unknown 
values  of  H0’,  and  hence  H0,  are  obtained  by  Gaussian 
elimination.  The  solution  H0  is  then  used  in  (4)  to 
obtain  the  electric  field  values  on  the  perimeter  of 
every  gridblock,  and  simple  averaging  is  used  to  obtain 
the  values  of  the  electric  fields  at  the  gridblock  centers. 

Typical  grids  consist  of  approximately  75  x  75 
gridblocks,  with  no  more  Uian  15  gridblocks  per 
wavelength.  Blocks  are  chosen  with  a  consideration  of 
the  radiation  wavelength,  the  radius  and  length  of  the 
antenna  being  modeled,  as  well  as  the  gap  size  used  to 
excite  the  antenna. 

Incorporation  of  stealth  boundary  conditions 

As  was  stated  earlier,  no  radiation  boundary  condition 
was  directly  applied  at  the  problem  domain's  outer 
surface.  Rather,  optimal  transmission  of  the  fields  into 
the  stealth  layer  was  sought.  The  tangential  electric 
field  at  the  outer  boundary  of  the  stealth  layer  was  set 
to  zero,  which  is  equivalent  to  bounding  the  stealth 
layer  with  a  perfect  conductor. 
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The  transmission  of  radiation  normally  incident  upon  a 
boundary  is  determined  by  the  wave  impedances  on 
both  sides  of  the  material  interface,  with  matched  wave 
impedances  providing  perfect  transmission.  To 
achieve  notm^  incidence  of  the  radiated  wave  at  the 
stealth  layer,  the  surface  of  the  stealth  layer  should  be  a 
sphere  centered  at  the  antenna,  with  a  radius  many 
times  the  wavelength.  In  practice  it  was  found 
sufficient  to  approximate  the  surface  of  the  sphere  by  a 
cylindrical  boundary,  as  shown  in  Fig.  2,  and  rely  upon 
multiple  reflections  from  the  stealth  layer  to  minimize 
the  disturbances  to  the  near  fields.  The  distance  from 
the  center  of  the  antenna  to  the  stealth  layer's  surface 
was  set  at  three  wavelengths.  A  chamfer  in  the  stealth 
layer,  as  indicated  by  the  dotted  line  in  the  upper  right 
hand  comer  of  Fig.  2,  was  initially  used  to  better 
approximate  a  spherical  boundary.  Only  a  small 
improvement  was  obtained,  however,  at  the  expense  of 
a  relatively  complex  design.  The  chamfer,  therefore, 
was  not  used  in  any  of  the  cases  considered  in  this 
paper. 

The  wave  impedance  and  propagation  constant  for  a 
plane  wave  propagating  in  a  general  lossy  medium  are 

V  e  ^}  (<T  +  jo)e) 


and  7  =  a  +  (8) 

where  f/stealth  =  ^problem  domain  would  maximize 
uansmission  for  normal  incidence  at  large  distances 
from  the  antenna  su^cture,  and  OCstealth  would  define 
the  attenuation  rate  in  the  absorbing  layer.  While  it  can 
be  seen  from  Fig.  2  that  waves  will  not  impinge  upon 
the  stealth  layer  normal  to  its  surface,  it  has  been  found 
that  choosing  the  electncal  properties  to  satisfy  the 
foregoing  provides  satisfactory  results,  as  exemplified 
by  the  case  studies  to  follow.  Given  an  excitation 
frequency,  p',  p",  a,  and  e  in  the  fictitious  stealth 
lay^  were  chosen  to  obtain  an  impedance  match,  and 
to  diminish  the  field  strength  to  insignificant  values  in 
the  allotted  stealth  layer  thickness.  When  choosing  the 
parameters  p',  p",  a,  and  e  for  the  stealth  layer,  a  great 
amount  of  flexibility  was  found  to  exist,  even  when 
specifying  both  the  stealth  layer  wave  impedance  and 
attenuation  per  wavelength.  The  stealth  layer  was 
typically  one  wavelength  (15  gridblocks)  thick, 
measured  in  the  stealth  material.  For  this  thickness, 
typically  substantially  less  than  1%  of  the  normally 
incident  power  was  reflected  from  the  surface  of  the 
stealth  layer,  and  the  absorbed  radiation  was  reduced  to 
approximately  1%  of  its  original  value,  before 
impinging  upon  the  perfectly  conducting  boundary 
which  is  assumed  to  bound  the  computational  grid. 
Thinner  stealth  layers  with  increased  attenuation  per 


wavelength,  or  thicker  stealth  layers  with  less 
attenuation  generally  also  provide  satisfactory  results. 
Our  own  choice  was  made  for  convenience  and  with  a 
view  to  not  significantly  increase  computauonal  times. 
For  example,  in  the  case  of  the  monopole  in  free  space 
0)  =  100  MHz,  p'  =  l.Opo,  p"  =  0.7po,  a  =  0.00389 
S/m,  and  e  =  1  .Ocq  in  the  stealth  layer.  This  matched 
the  free  space  wave  impedance  and  provided  98.8% 
attenuation  per  wavelength. 

Simulation  results 

This  section  presents  some  typical  results  of  the  thin 
wire  antenna  configurations  modeled.  Simulation 
results  are  compared  to  both  experimental  and 
theoretical  data.  HERZ  was  implemented  in 
FORTRAN  using  Macintosh  Programmer's  Workshop 
(MPW)  version  3.0  on  a  Macintosh  II  personal 
computer,  with  an  RP88  coprocessor  to  increase 
performance,  which  resulted  in  execution  times  of 
approximately  six  minutes. 

Three  representative  antenna  configurations  were 
considered:  (1)  A  bare  coaxially  fed  monopole  over  an 
infinite  conducting  plane  operating  in  a  lossless 
medium,  (2)  a  bare  center  fed  dipole  operating  in  a 
lossy  medium,  and  (3)  a  dielectrically  coated  monopole 
operating  in  air.  Fig.  4  displays  the  general 
configurations  modelled  in  each  case. 

To  validate  HERZ,  the  surface  current  distribution  and 
input  admittance  for  each  antenna  configuration  were 
compared  to  available  theoretical  or  measured  values. 
The  boundary  condition  for  the  tangential  magnetic 
field  at  the  antenna's  surface,  nxH  =  J^,  was  used  to 
determine  the  surface  current  I,  directly  from  H0  = 
IfZitr.  The  driving  point  admittance  was  obtained  by 
dividing  the  current  on  the  antenna  nearest  the  driving 
point  by  the  voltage  across  the  gap  at  that  point, 
calculated  from  jE-dl,  where  dl  is  perpendicular  or 
parallel  to  the  antenna  axis,  depending  on  whether  the 
gap  is  coaxial  or  that  of  a  center  fed  dipole.  Fig.  4 
displays  the  driving  points  for  all  configurations. 
HERZ  calculates  the  fields  at  the  gridblock  centers,  so 
to  find  H0  on  the  antenna's  surface,  a  third-order 
polynomial  was  fitted  to  the  field  values  in  the 
neighboring  blocks,  extrapolating  radially  to  the 
surface.  The  gridblock  configuration  in  the  excitation 
gaps  will  be  given  for  each  case. 

Fig.  5(a)  displays  the  driving  point  admittance  versus 
antenna  length  for  the  monopole  in  free  space.  This 
thin  antenna,  with  aJX  =  0.0064,  and  b/a  =  1.189,  was 
excited  at  100  MHz.  The  theoretical  curve,  as 
determined  by  King  (1971,  p.  11),  was  matched  very 
well  by  the  output  from  HERZ,  for  all  antenna  lengths. 
Fig.  5(b)  displays  the  normalized  current  distribution 
along  a  similar  monopole  in  free  space,  with  h/X.  =  0.5 
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and  aJX  =  0.0254,  also  excited  at  100  MHz.  Here, 
HERZ  agrees  closely  with  the  transmission  line  model 
by  King  (1971,  p.  18).  For  both  cases,  the  excitation 
was  specified  at  the  ground  plane  level,  assuming  a  1/r 
coaxial  voltage  distribution  across  the  gap,  with  only  2 
gridblocks  across  this  region.  The  formulas  used  in 
King's  theory  (1971,  p.  9)  rest  on  the  assumption  that  a 
TEM  mode  exists  at  the  junction  of  the  coaxial  feed 
and  the  antenna,  justifying  the  1/r  excitation  used  in 
HERZ. 

Fig.  6(a)  shows  the  excellent  agreement  between 
HERZ's  calculated  admittance  and  Scott's  measured 
values  and  King's  theorei  oal  results  (1981  p.  170)  for  a 
dipole  radiating  in  a  lossy  medium,  with  a  conductivity 
of  5.35x10'^  S/m,  and  relative  dielectric  constant  of 
6.0.  For  this  trial,  a/k  =  0.00265,  and  a/b  =0.07,  at  1 14 
MHz,  where  X,  a,  and  P  are  the  wavelength, 
attenuation  constant,  and  propagation  constant  of  the 
radiation  in  the  lossy  medium  (Ramo  and  Whinnery 
1967).  Examination  of  Fig.  6(a)  shows  that  HERZ 
provides  closer  agreement  with  the  measured  values 
than  do  King's  theoretical  results. 

King  assumed  a  delta  function  generator  in  an 
infinitesimal  gap  in  formulating  his  theoretical  curve, 
but  unfortunately,  Scott's  gap  size  is  not  specified  in  the 
forementioned  reference.  HERZ  utilized  a  gap  size 
(1.0  mm)  of  the  same  order  as  the  antenna  radius  (2.65 
mm)  for  this  trial,  but  several  other  trials  indicate  that 
the  gap  size  is  not  critical,  with  results  remaining 
nearly  constant  for  all  gap  sizes  less  than  the  antenna 
radius.  The  critical  factor  in  modeling  the  gap  was  the 
number  of  blocks  employed,  with  four  radial  and  three 
axial  blocks  sufficiently  representing  this  sensitive 
region  for  the  trial  represented,  as  the  good  results 
indicated.  As  the  number  of  gridblocks  in  this  region 
was  decreased,  the  admittance  results  deviated  from 
Scott's  measured  values,  but  increasing  the  number  of 
blocks  increased  the  simulation  run  lime,  while  having 
insignificant  effects  on  the  simulation  results.  The 
normalized  current  distribution  for  this  antenna 
configuration  was  also  determined  using  HERZ,  for  an 
antenna  of  length  ph  =  0.315.  Fig.  6(b)  shows  the 
simulation  results  versus  both  Scott's  measured  values 
and  King's  theoretical  distribution  (1981,  p.  165). 

Experimental  results  for  dielecmcally  coated  antennas 
by  Lamensdorf  (1967)  were  used  to  verify  HERZ  in 
this  capacity.  Figures  7(a)  and  7(b)  display 
Lamensdorfs  measured  input  admittance  for  a 
monopole  operating  in  free  space  at  600  MHz  with  2a 
=  6.35  mm,  2b  =  19.05  mm,  and  D/2a  =  3.74,  for  two 
different  dielectric  coalings,  Cr  =  9.0  and  tj  =  15.0. 
Both  dielectric  coatings  were  modeled  with 
conductivities  of  0.001  S/m.  This  value  was  chosen  to 
best  represent  the  non-specific  value  of  a  <  0  001  S/m 
provided  in  the  Lamensdorf  publication.  HERZ’s 


simulation  results  replicate  the  shape  of  Lamensdorfs 
measured  distribution,  providing  a  reasonable  overall 
match  in  both  cases.  Unlike  the  monopole  previously 
modeled,  in  which  a  TEM  coaxial  field  distribution 
across  the  gap  could  be  assumed,  the  coaxial  line  was 
modeled  to  best  represent  the  actual  fields  across  the 
gap  at  the  ground  plane  level.  The  number  of  radial 
gridblocks  used  to  model  the  coaxial  line  feeding  the 
antenna  depended  on  the  ratio  of  the  outer  to  inner 
coaxial  radii,  b/a.  The  previous  monopole  modeled  had 
a  b/a  value  of  1.189  so  two  grid  blocks  to  represent  the 
gap  were  sufficient.  In  the  present  example,  b/a  =  3.0, 
and  trial  and  error  indicated  that  six  radial  gridblocks 
were  optimal  in  modeling  the  gap.  TEM  excitation  was 
applied  at  one  eighth  of  a  wavelength  of  the  coaxial 
line  below  the  ground  plane.  This  length  of  line  was 
modeled  in  ten  axial  gridblocks.  HERZ's  current 
distribution  on  a  similar  antenna  configuration  with  2a 
=  6.35  mm,  2b  =  19.05  mm  and  a  =  0.0032  S/m,  but 
with  D/2a  =  8,  is  displayed  in  Fig.  7(c)  along  with 
Lamensdorfs  measur^  values.  The  agreement  is 
reasonable,  but  not  as  good  as  in  previous  cases.  This 
may  reflect  the  difficulty  of  modeling  a  larger  coaxial 
gap  than  previously  considered,  or  the  inhomogeneity 
associated  with  the  long  dielectric  coaling. 

Additional  antenna  configurations,  radiating  into 
materials  ranging  from  air  to  highly  lossy,  were  tested 
along  with  those  cases  considered  in  this  paper. 
Similar  agreement  with  theoretical  and  measured  data 
was  obtained.  This  agreement  was  generally 
comparable  to  or  better  than  that  achieved  earlier  with  a 
finite  element  program  utilizing  radiation  boundary 
conditions  (Sumbar  ei  a/.  1990). 

Conclusions 

The  use  of  stealth  boundary  conditions  in  the  finite 
difference  simulation  of  simple  antenna  problems 
appears  to  be  very  effective.  The  method  is  physically 
very  transparent  and  the  simplicity  of  the  finite 
difference  grid  suucture  permits  easy  realization  of  the 
antenna  configurations.  The  methodology  was  tested 
for  a  monopole  in  free  space,  a  dipole  in  a  lossy 
medium,  and  a  dielectrically  coated  monopole  in  free 
space.  In  all  cases,  results  agreed  well  with  published 
theoretical  and  experimentally  measured  values. 
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Figure  1 

The  general  antenna  configuration  and 
coordinate  system  used  in  HERZ. 
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Figure  2 

The  outline  of  the  finite  difference  grid  system  used  to  model  the  problem  domain  and  stealth 
layer.  All  boundary  conditions  are  specified  on  the  outermost  boundary,  as  indicated. 
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Figure  3 

The  finite  difference  grid  stucture  used  in  HERZ.  All  magnetic  field  components  are 
azimuthally  directed  and  represented  at  the  block  centers,  while  the  electric  field 
components  are  represented  on  the  block  boundaries. 
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Figure  S<b)  Monopole  in  air 

The  normalized  current  distribution 
for  a  monopole  on  a  conductive  plane 
radiating  into  a  lossless  medium. 
Theoretical  data  published  by  King 
(1971,  p.ll)  is  used  to  validate  HERZ's 
calculated  values. 
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Figure  6(a)  Dipole  in  a  lossy  medium 

The  driving  point  admittance  versus 
antenna  length  for  a  dipole  operating 
in  a  lossy  medium  as  measured  by  Scott, 
theoretically  determined  by  King  (1981, 
p.  170),  and  calculated  with  HERZ. 
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antenna  admittance  (mS)  position  on  antenna  (z/h) 
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Figure  6(b)  Dipole  in  a  lossy  medium 

HERZ's  normalized  current  distribution 
for  a  dipole  radiating  into  a  lossy 
medium  is  compared  to  theoretical 
data  by  King  (1981,  p.l6S)  and  the 
measured  values  by  Soon. 
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Figure  7(a)  Dielectrically  coated 
monopole  in  air 

The  driving  point  admittance  versus 
antenna  length  for  a  dielectrically 
coaled  monopole  in  air  as  measured 
by  Lamesndorf  (1967)  and  calcualted 
with  HERZ. 
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Figure  7(b)  Dielectrically  coated 
monopole  in  air 

The  driving  point  admittance  versus 
antenna  length  for  a  dielectrically 
coated  monopole  in  air  as  measured 
by  Lamensdorf  (1967)  and  calculated 
by  HERZ. 
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Figure  7(c)  Dielectrically  coated 
monopole  in  air 

The  normalized  current  distribution  for 
a  dielectrically  coated  monopole  on  a 
conducting  plane  operating  in  air.  HERZ's 
calculated  values  are  compared  to 
Lamensdorf  s  measured  values  (1967). 
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ABSTRACT 

The  two  wire  rhombic  illuminator  (RI)  is  a  rhombic  wire  antenna  located  over  a  ground  plane  which 
is  used  to  generate  near  un^orm  electromagnetic  fields  over  a  specified  working  volume  located  below 
the  antenna.  The  RI  can  be  modeled  as  a  single  two-wire  transmission  line  at  low  frequencies  but  the 
transmission  line  model  becomes  inadequate  at  higher  frequencies  due  to  the  enhanced  radiation 
properties  of  the  antenna.  The  performance  of  the  RI  including  radiation  effects  is  evaluated  numerically 
using  the  Numerical  Electromagnetics  Code  (NEC-2)  and  compared  to  that  of  the  transmission  line  model. 
The  antenna  in  the  presence  of  a  perfectly  conducting  ground  plane  is  analyzed  at  frequencies  ranging 
from  I  MHz  to  500MHz.  Various  segmentation  schemes  are  utilized  in  order  to  accurately  predict  the 
performance  of  the  RI  at  frequencies  where  the  dimensions  of  the  structure  are  very  large  in  terms  of 
wavelength. 


INTRODUCTION 

The  rhombic  transmission  line  configuration  has  received  much  attention  in  recent  years  as  an 
illuminator  for  surveillance  testing  of  electromagnetic  shielding.  The  first  detailed  study  of  such  an 
illuminator  was  performed  by  Baum  [1]  and  Shen  and  King  [2,3,4].  The  uniformity  of  the  working 
volume  fields  for  the  two-wire  rhombic  illuminator  were  andyzed  by  Zuffada  and  Engheta  [S]  using  a 
two-wire  transmission  line  operating  in  the  TEM  mode.  As  noted  in  [5],  the  accuracy  of  the  transmission 
line  model  diminishes  as  the  frequency  is  increased  due  to  the  appearance  of  longitudinal  modes. 

In  this  paper,  the  validity  of  the  RI  transmission  line  model  is  investigated  using  a  numerical  solution 
technique.  An  illustration  of  the  RI  over  a  perfectly  conducting  ground  plane  is  shown  in  Figure  1 .  The 
actual  RI  dimensions  and  working  volume  location  are  shown  in  Figure  2.  All  wires  have  a  diameter  of 
3.175  mm.  The  anterma  current  and  the  individual  components  of  the  working  volume  electric  field  are 
computed  over  a  range  of  1  MHz  to  500 MHz.  Sever^  aspects  of  the  RI  performance  are  noted  with 
regard  to  the  behavior  of  the  antenna  currents  and  working  volume  fields  with  frequency. 

The  RI  may  be  driven  in  the  common-mode  (push-push)  configuration  with  respect  to  the  ground 
plane.  The  four  wires  of  the  RI  are  denoted  as  wire  #1  and  wire  #3  (which  make  up  the  launch  region) 
which  are  connected  with  wire  #2  and  wire  #4,  respectively,  (which  make  up  the  termination  region). 
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Figure  2.  Rhombic  illuminator  and  working  volume  dimensions. 
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The  common  mode  excitation  of  the  RI  yields  an  electric  field  which  is  predominantly  vertical  in  the 
given  working  volume. 


TRANSMISSION  LINE  MODEL  OF  THE  RHOMBIC  ILLUMINATOR 

As  previously  mentioned,  the  RI  may  be  modeled  by  a  two-wire  traitsmission  line  over  a  perfectly 
conducting  ground  plane  at  low  frequencies.  The  TEM  fields  of  the  two-wire  line  under  common-mode 
excitation  are  obtained  by  considering  the  static  field  of  equivalent  line  charges  over  a  perfectly  conducting 
ground  plane.  For  the  given  woilung  volume  dimensions,  the  values  of  the  wire  separation  and  wire 
height  above  ground  may  be  optimized  with  regard  to  field  uniformity  [5], 

The  larger  vertical  component  of  the  common-mode  electric  field  is  designated  as  the  principal 
component  while  the  smaller  horizontal  electric  field  component  is  referred  to  as  the  nonprincipal  field 
component.  The  principal  and  nonprincipal  fields  (£,  and  E,,  respectively)  for  the  common-mode 
excitation  [5,8,9]  are  given  by 
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where  q  is  the  line  charge  density,  is  the  permittivity  of  vacuum,  la  is  the  distance  between  the  wire 
centers  and  b  is  the  height  of  the  wires  above  the  ground  plane.  The  value  of  the  line  charge  density  in 
(1)  and  (2)  is  related  to  the  potential  difference  between  each  of  the  wires  and  ground  (V)  and  the 
characteristic  impedance  of  the  common-mode  two-wire  line  shown  in  Figure  3  as  given  by 


_?_=30_L 
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The  characteristic  impedance  of  the  two-wire  line  (Z^*)  under  common-mode  excitation  is 
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where  Tjo  is  the  intrinsic  impedance  of  vacuum  and  r  is  the  radius  of  the  wires. 


V  N  b 
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Figure  3.  Common-mode  two  wire  line  over  ground. 
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The  principal  and  nonprincipal  transverse  electric  fields  in  the  working  volume  of  the  RJ  are 
symmetric  about  the  y-z  plane  given  common-mode  excitation.  Therefore,  the  field  behavior  over  the 
entire  transverse  cross-section  of  the  working  volume  given  by  [-22.5m<x<22.5m]  and  [0<z<  15m]  is 
evident  given  the  fields  over  the  cross-section  defined  by  [0<x<22.5m]  and  [0<z<  15m].  The  transverse 
cross-section  located  at  the  center  of  the  working  volume  (y=125m)  is  designated  as  and  is  shown  in 
Figure  4.  For  purposes  of  comparison,  the  individual  electric  field  components  are  plotted  over  .  The 
dimensions  of  the  RI  at>’=  125m  are  2a=25.48m  (wire  separation)  and  /)=20.70m  (wire  height  above  the 
ground  plane). 


Figure  4.  Working  volume  cross-sectional  surface  (S^)  over  which  fields  are  plotted. 

The  transverse  electric  field  components  of  a  similarly  spaced  two-wire  transmission  line  over  a 
ground  plane  (2a=25.48m,  A=20.70m,  r=  1.5875mm)  with  V=  1  V  are  shown  in  Figures  5  and  6.  The 
principal  electric  field  component  of  the  transmission  line  model  shown  in  Figure  5  contains  a  broad  peak 
along  the  upper  edge  of  the  working  volume  (z=15m)  which  corresponds  to  points  directly  below  the 
transmission  line  conductor.  The  principal  electric  field  is  relatively  uniform  over  the  remainder  of  the 
working  volume.  The  nonprincipal  electric  field  component  shown  in  Figure  6  is  zero-valued  along  the 
lower  edge  of  the  working  volume  (tangent  to  the  ground  plane)  and  increases  in  magnitude  as  the 
observation  point  is  moved  away  from  the  ground  plane.  The  nonprincipal  electric  field  contains  a  null 
along  the  upper  edge  of  the  working  volume  which  corresponds  to  the  point  where  the  horizontal 
components  due  to  both  wires  and  their  images  cancel  one  another.  The  x-coordinate  of  this  null  lies 
between  the  transmission  line  center  point  at  x=0  and  the  coordinate  of  the  transmission  line  conductor 
at  x=a.  The  transverse  fields  of  the  transmission  line  model  are  compared  to  the  fields  generated  by  the 
corresponding  NEC  thin- wire  model  of  the  RI  described  in  the  following  section. 

NUMERICAL  MODEL  OF  THE  RHOMBIC  ILLUMINATOR 

The  number  of  segments  required  to  model  each  of  the  four  wires  of  the  RI  using  NEC-2  is  based  on 
the  frequency  of  operation.  In  general,  the  segment  lengths  should  be  chosen  in  the  range  of  0.00  IX  to 
O.IX  to  maintain  numerical  stability  [6].  This  segment  length  requirement  makes  the  analysis  of 
electrically  long  thin-wire  structures  (such  as  the  RI  at  high  frequencies)  computationally  prohibitive. 
However,  the  segment  lengths  can  be  exponentially  graded  over  the  entire  wire  ("entire  grading")  or  a 
portion  of  each  wire  ("partial  grading")  [7,8,9]  to  minimize  the  number  of  segments  required  for  an 
accurate  solution  of  the  RI  geometry.  These  segmentation  schemes  allow  for  electrically  short  segments 
in  critical  regions  (source,  termination,  wire  bends)  while  allowing  for  longer  segments  over  the  remainder 
of  the  structure. 
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Figure  5.  Principal  transverse  electric  field  Figure  6.  Nonprincipal  transverse  electric  field 
magnitude  on  (transmission  line  model).  magnitude  on  (transmission  line  (model). 


The  segments  on  the  launching  region  wires  (#1  and  #3)  are  numbered  from  the  ground  plane  to  the 
apex  of  the  RI.  That  is,  the  first  segment  on  each  launch  wire  is  connected  to  the  ground  plane  while  the 
last  segment  is  connected  to  the  termination  wire  at  the  apex.  The  segments  on  the  termination  region 
wires  (#2  and  #4)  are  numbered  from  the  apex  to  the  ground  plane.  The  two  common-mode  voltage 
sources  (applied  electric-field  sources),  and  Vs:,  are  placed  on  the  second  segment  of  wire  #1  and  wire 
#3,  respectively.  The  magnitude  of  both  voltage  sources  is  one  volt  with  a  phase  of  zero  degrees.  The 
terminations  placed  on  the  last  segments  of  wire  #2  and  wire  #4  are  lumped  resistive  loads  of  630 Q.  The 
630 n  lumped  resistive  value  corresponds  to  twice  the  average  common-mode  characteristic  impedance 
along  the  length  of  the  Rl  as  given  in  Equation  (4). 

NEC  allows  for  the  utilization  of  structural  symmetry  to  simplify  the  geometry  input  as  well  as  the 
solution  process.  The  RJ  configuration  shown  in  Figure  1  is  symmetric  about  the  y-z  plane.  Thus,  using 
symmetiy,  only  wire  #1  and  wire  #2  are  needed  to  describe  the  entire  RI  structure.  The  voltage  sources 
are  not  reflected  using  symmetry  so  that  both  sources  of  the  RI  must  be  defined  explicitly.  NEC  does, 
however,  allow  for  the  reflection  of  loaded  segments  so  that  the  resistively-loaded  segment  on  wire  #2  is 
automatically  reflected  to  wire  #4  using  symmetry. 

The  previously  defined  thin-wire  model  of  the  RJ  is  analyzed  using  discrete  NEC  models  over  a 
frequency  range  of  1  MHz  to  500  MHz.  The  NEC  models  developed  for  the  RI  are  highly  sensitive  to  the 
frequency  of  operation  because  of  the  critical  nature  of  the  segment  length  with  regard  to  wavelength. 
For  this  reason,  the  RI  model  is  analyzed  at  five  distinct  frequencies:  I  MHz,  10  MHz,  100  MHz,  300  MHz 
and  500MHz.  The  RI  current  and  electric  field  components  over  are  computed  at  each  frequency. 
The  total  number  of  unknowns  iN,„^i)  is  equal  to  the  total  number  of  segments  required  to  model  a 
symmetric  cell  (wire  #1  plus  wire  #2).  Within  NEC,  the  maximum  array  dimension  with  regard  to  the 
number  of  segments  must  be  2N,„^i.  The  total  length  of  the  RI  symmetric  cell  is  approximately  265  m. 
The  total  electrical  lengths  of  the  RI  symmetric  cell  at  the  various  frequencies  of  interest  are  listed  in 
Table  1. 

The  RI  at  1  MHz  is  easily  modeled  with  "optimum”  length  segments  of  0.0  IX.  The  length  of  0.01  X 
lies  one  order  of  magnitude  above  the  recommended  minimum  segment  length  and  one  order  of  magnitude 
below  the  maximum  recommended  segment  length.  Wire  #1  and  wire  #2  are  divided  into  equal  segments 
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of  length  approximately  O.OIX.  which  yields  54  segments  on  wire  #1  and  35  segments  on  wire  #2  for  a 
total  of  89  segments  in  a  symmetric  cell. 

The  total  length  of  wire  #1  plus  wire  #2  at  lOMHz  is  approximately  8.8X.  Thus,  approximately  880 
segments  would  be  required  to  model  the  RJ  with  O.OIX.  segments.  Yet,  the  use  of  this  many  segments 
yields  a  solution  which  is  unnecessarily  lengthy.  The  accuracy  of  the  NEC  model  for  long  straight  wires 
is  not  appreciably  affected  by  using  segments  which  are  longer  than  0.01  X.  away  from  the  critical  regions 
of  the  structure  [6].  Thus,  by  exponentially  grading  the  segments  along  the  RI  wires,  the  total  number 
of  segments  required  can  be  reduced  significantly.  The  segments  near  the  source  must  be  short  and  of 
comparable  length  to  yield  accurate  representations  of  the  applied  electric  field  sources.  Therefore,  the 
first  segment  on  wire  #1  is  chosen  to  be  O.OIX.  in  length.  The  partial  grading  technique  is  then  used  with 
a  grading  factor  of  1.1  (10%  increase  in  length  from  segment  to  segment)  until  the  segments  reach  a 
maximum  length  of  A„,^,=0.1X..  The  same  segment  grading  is  used  at  the  apex  end  of  wire  #1.  The 
remaining  midsection  of  wire  #1  is  modeled  with  uniform  length  segments  of  length  A^.  The  same 
segmentation  scheme  is  applied  to  the  termination  wire  (wire  #2)  resulting  in  83  segments  on  wire  #1  and 
65  segments  on  wire  #2  for  a  total  of  148  segments. 

At  100  MHz,  the  total  electrical  length  of  the  symmetric  cell  is  approximately  88X,.  In  order  to 
significantly  reduce  the  number  of  segments  required  at  this  frequency,  the  partial  grading  scheme  must 
be  employed.  As  in  the  10  MHz  case,  the  shortest  segments  on  wire  #1  and  wire  #2  are  defined  with 
lengths  of  O.OIX.  and  a  grading  factor  of  1 . 1  is  used  for  the  weighted  portions  of  the  wires.  The  maximum 
segment  length  is  chosen  to  be  A,^=0.2X..  This  segmentation  results  in  310  segments  on  wire  #1  and  218 
segments  on  wire  #2  which  yields  A,„„,=528. 

The  total  electrical  length  of  the  symmetric  cell  at  300 MHz  is  approximately  264X..  The  same 
segmentation  scheme  as  outlined  for  the  lOMHz  and  lOOMHz  cases  is  used  for  the  300 MHz  case.  The 
grading  factor  is  again  cnosen  to  be  1.1  while  the  maximum  segment  length  is  chosen  as  0.3X..  The 
resulting  segment  totals  on  wire  #1  and  wire  #2  are  585  and  400,  respectively.  Thus,  the  total  number 
of  segments  in  the  symmetric  cell  is  985. 

The  previously  defined  segmentation  scheme  is  applied  to  the  RI  at  500  MHz.  The  electrical  length 
of  the  RI  symmetric  cell  at  500  MHz  is  in  excess  of  440X..  A  grading  factor  of  1.1  and  a  maximum 
segment  length  of  0.4X,  are  utilized.  This  model  of  the  RI  yields  724  segments  on  wire  #1  and  493 
segments  on  wire  #2.  Thus,  the  entire  symmetric  cell  requires  a  total  of  1217  segments. 

The  characteristics  of  the  five  NEC  models  are  summarized  in  Table  1.  The  use  of  NEC  models 
which  contain  segment  lengths  in  excess  of  O.IX.  requires  special  consideration  with  regard  to  the 
interaction  approximation  distance  within  NEC.  When  segments  are  separated  by  more  than  a  prescribed 
distance  (known  as  the  interaction  approximation  distance),  the  source  segment  is  modeled  as  an 
infinitesimal  current  element  at  the  segment  center.  This  approximation  is  accurate  when  using  short 
segments.  However,  longer  segments  may  have  significant  variation  of  current  magnitude  and  phase  along 
the  length  of  the  segment  which  makes  this  approximation  inaccurate.  Thus,  when  using  longer  segments, 
the  interaction  approximation  must  not  be  utilized.  The  interaction  approximation  distance  may  be  defined 
by  the  user  but  a  default  value  of  1  m  is  assumed  in  NEC.  In  the  results  given  here,  the  interaction 
approximation  distance  is  defined  large  enough  that  the  interaction  approximation  is  not  implemented. 

NUMERICAL  RESULTS  FOR  THE  THIN-WIRE  RI  MODEL 

The  magnitude  of  the  resulting  current  along  wire  #1  and  wire  #2  of  the  RI  is  shown  in  Figure  7.  The 
current  magnitude  at  1  MHz  exhibits  a  standing  wave  pattern  that  is  very  slightly  damped  along  the  length 
of  the  two  wires.  The  principal  transverse,  nonprincipal  transverse  and  longitudinal  components  of  the 
electric  field  over  are  shown  in  Figures  8,  9  and  10,  respectively.  The  transverse  electric  field  plots 
may  be  compared  with  those  of  the  two-wire  transmission  line  (Figures  5  and  6).  As  expected,  the  general 
distribution  of  transverse  electric  field  components  are  nearly  identical  to  those  of  the  two-wire 
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(MHz) 

Total  length 
(wire  #1  +  wire  #2) 

Segment  grading 
scheme 

^max 

1 

0.8822X 

none 

0.0  IX 

89 

10 

8.822X 

partial 

O.IX 

148 

100 

88.22A. 

partial 

0.2X 

528 

300 

264.7X 

partial 

0.3X 

985 

500 

441. IX 

partial 

0.4X 

1217 

Table  1.  Descriptions  of  the  R1  NEC  models. 

transmission  line  TEM  mode.  The  magnitudes  of  the  transverse  electric  field  components  for  the  RI 
model  are  slightly  higher  than  the  corresponding  transmission  line  values.  The  magnitude  of  the 
longitudinal  electric  field  component  is  small  in  comparison  with  the  transverse  components. 

The  magnitude  of  the  current  along  wires  #1  and  #2  of  the  RI  at  10  MHz  is  plotted  in  Figure  1 1 .  The 
corresponding  electric  field  components  are  plotted  over  and  appear  in  Figures  12,  13  and  14.  The 
magnitude  of  the  RI  source  current  at  10  MHz  exhibits  a  significant  increase  when  compared  to  the  1  MHz 
value.  Yet,  the  RI  current  at  lOMHz  exhibits  a  slightly  higher  rate  of  decay  than  the  1  MHz  distribution. 
The  general  shapes  of  the  electric  field  components,  both  transverse  and  longitudinal,  remain  unchanged 
from  the  IMHz  values.  The  magnitudes  of  all  three  electric  field  components  at  10  MHz  are  slightly 
larger  than  the  corresponding  1  MHz  values. 

The  100 MHz  RI  current  magnitude  is  shown  in  Figure  15.  The  resulting  100  MHz  working  volume 
electric  field  components  are  given  in  Figures  16,  17  and  18.  An  increase  in  the  magnitude  of  the  source 
current  from  10  MHz  to  100  MHz  is  observed.  However,  the  damping  factor  of  the  100  MHz  current  is 
larger  than  that  of  the  lOMiIz  current  yielding  only  slightly  higher  values  of  current  magnitude  in  the 
vicinity  of  the  working  volume  (near  the  wire  apex).  The  magnitudes  of  both  the  transverse  and 
longitudinal  components  of  the  working  volume  electric  field  exhibit  a  slight  increase  over  the 
corresponding  values  at  10  MHz.  The  variation  in  the  working  volume  field  components  with  respect  to 
wavelength  becomes  apparent  at  lOOMHz  as  shown  in  Figures  16,  17  and  18. 

The  current  magnitude  along  the  RI  wires  at  300 MHz  is  plotted  in  Figure  19,  while  the  corresponding 
electric  field  components  over  are  plotted  in  Figures  20,  21  and  22.  A  general  trend  is  evident  in  the 
variation  of  the  RI  current  and  working  volume  field  values  as  a  function  of  frequency.  As  the  frequency 
of  operation  is  increased,  the  source  current  magnitude  is  increased.  The  standing  wave  pattern  along  the 
RI  wires  stays  relatively  constant  with  frequency.  The  resulting  working  volume  electric  field  components 
exhibit  a  small  increase  in  magnitude  as  the  frequency  is  increased. 

The  magnitude  of  the  current  on  wires  #1  and  #2  of  the  RI  at  500MHz  is  shown  in  Figure  23.  The 
electric  field  components  in  the  working  volume  of  the  RI  at  500  MHz  are  shown  in  Figures  24,  25  and 
26.  The  current  magnitude  at  the  source  is  slightly  larger  than  that  of  the  300 MHz  case  but  the  current 
distribution  away  from  the  source  is  quite  similar  to  that  of  the  300  MHz  case.  The  standing  wave 
amplitude  is  slightly  higher  than  that  of  the  300  MHz  case.  The  working  volume  field  components  exhibit 
a  small  increase  in  amplitude  over  the  300  MHz  case  as  expected. 

In  order  to  test  the  convergence  of  the  previous  NEC  solutions  when  using  relatively  long  segments, 
the  300  MHz  and  500  MHz  cases  are  modeled  with  slightly  shorter  segments  and  the  resulting  currents  are 
compared  with  the  currents  obtained  with  the  initial  models.  These  comparisons  are  shown  in  Figures  27 
and  28.  For  both  the  300 MHz  and  500 MHz  cases,  the  current  distributions  are  quite  similar  with  only 
minor  variations  in  the  standing  wave  amplitude. 
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Figure  11.  RI  current  magnitude  on 
wires  #1  and  #2  at  10  MHz. 


Figure  12.  Principal  transverse  electric 
field  magnitude  on  5^  at  10  MHz. 


Figure  13.  Nonprincipal  transverse  electric 
field  magnitude  on  at  10  MHz. 


Figure  14.  Longitudinal  electric  field 
magnitude  on  at  10  MHz. 
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Figure  27.  Comparison  of  RI  current  magnitude  at  300MHz  with  A^=0.3X  and  A^=0.2X 
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Figure  28.  Comparison  of  RI  current  magnitude  at  500 MHz  with  A^=0.4A,  and  A^=0.3X 

ANALYSIS  OF  RESULTS  AND  CONCLUSIONS 

The  rhombic  illuminator  configuration  has  been  analyzed  using  a  numerical  solution  technique  over 
the  frequency  range  of  IMHz  to  500  MHz.  The  numerical  solution  is  obtained  using  NEC  with  five 
discrete  models  at  1,  10,  100,  300  and  500 MHz.  The  antenna  current  and  the  working  volume  electric 
field  components  were  computed  and  compared  to  the  fields  obtained  using  the  two-wire  transmission  line 
TEM  model.  The  current  along  the  wires  of  the  illuminator  was  found  to  display  several  basic  physical 
characteristics.  The  current  exhibited  a  damped  standing  wave  pattern  which  remains  relatively  constant 
over  the  frequency  range  of  interest.  Reflections  from  the  wire  bend  become  evident  at  frequencies  of 
100  MHz  and  above.  The  current  magnitude  in  the  source  region  increases  rapidly  from  1  to  100  MHz 
and  at  a  slower  rate  at  the  higher  frequencies.  Overall,  the  average  current  magnitude  along  the  RI  wires 
is  relatively  constant  over  the  given  frequency  range. 
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The  electric  field  components  in  the  working  volume  also  exhibit  several  distinct  characteristics.  The 
magnitudes  of  the  working  volume  fields  increase  with  frequency  since  the  peak  currents  on  the  antenna 
also  increase  with  frequency.  The  general  shapes  of  the  field  component  distributions  depend  on  the 
distribution  of  current  on  the  illuminator  wires.  Since  the  current  magnitude  distribution  is  relatively 
constant  over  the  frequency  range  of  interest,  the  field  component  distributions  also  remain  relatively 
constant.  The  magnitude  of  the  longitudinal  electric  field  component  remains  quite  small  in  comparison 
to  the  principal  transverse  component  indicating  that  no  significant  longitudinal  modes  are  excited  at 
frequencies  of  500MHz  and  below.  As  shown  in  Figures  8,  12,  16,  20  and  24,  the  R1  performance  is 
satisfactory  up  to  500  MHz  with  regard  to  the  uniformity  of  the  principal  transverse  field  in  the  working 
volume.  The  transmission  line  model  yields  an  accurate  prediction  of  the  general  shape  of  the  working 
volume  transverse  field  component  distributions  but  does  not  accurately  predict  the  magnitudes  of  these 
fields  [10].  At  500  MHz,  the  transverse  fields  are  roughly  three  times  larger  than  those  predicted  by  the 
theoretical  model.  It  may  be  possible  to  reduce  the  standing  wave  amplitude  along  the  structure  at  the 
various  frequencies  by  improved  load  matching  or  through  the  use  of  distributed  loading  along  the  wires. 
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PUBLICATION  CRITERIA 

Each  paper  is  required  to  manifest  some  relation  to  applied  computational  electromagnetics. 
Papers  mav  address  general  issues  in  applied  computational  electromagneUcs.  or  they  mav  focus  on 
specific  applications,  techniques,  codes,  or  computational  issues.  While  the  following  list  is  not 
exhaustive,  each  paper  will  generadly  relate  to  at  least  one  of  these  areas: 

1.  Code  validation.  This  is  done  using  internal  checks  or  experimental,  analytical  or  other 
computational  data.  Measured  data  of  potential  utility  to  code  validation  efforts  will  also  be  considered 
for  publication. 

2 .  Code  performance  analysis .  This  usually  involves  identification  of  numerical  accuracy  or  other 
limitations,  solution  convergence,  numerical  and  physical  modeling  error,  and  parameter  tradeoffs. 
However,  it  is  also  permissible  to  address  Issues  such  as  ease-of-use,  set-up  time,  run  time,  special 
outputs,  or  other  special  features. 

3 .  Computational  studies  of  basic  physics.  This  Involves  using  a  code,  algorithm,  or  computational 
technique  to  simulate  reality  in  such  a  way  that  better  or  new  physical  Insight  or  understanding  is 
achieved. 

4.  New  computational  techniques,  or  new  applications  for  existing  computational  techniques  or 

codes. 

5.  “Tricks  of  the  trade"  In  selecting  and  applying  codes  and  techniques. 

6.  New  codes,  algorithms,  code  enhancement,  and  code  fixes.  This  category  is  self-explanatory 
but  includes  significant  changes  to  existing  codes,  such  as  applicability  extensions,  algorithm  optimiza¬ 
tion,  problem  correction,  limitation  removal,  or  other  performance  improvement.  NOTE:  CODE  (OR 
ALGORITHM)  CAPABILITY  DESCRIPTIONS  ARE  NOT  ACCEPTABLE.  UNLESS  THEY  CONTAIN  SUFH- 
CIENT  TECHNICAL  MATERIAL  TO  JUSTIFY  CONSIDERATION. 

7.  Code  Input/output  Issues.  This  normalfy involves  innovations  in  input  (such  as  input  geometry 
standardization,  automatic  mesh  generation,  or  computer-aided  design)  or  in  output  (whether  it  be 
tabular,  graphical,  statistical,  Fourier- transformed,  or  otherwise  signal-processed).  Material  dealing  with 
Input/output  data  base  management,  output  interpretation,  or  other  input/ output  Issues  will  also  be 
considered  for  publication. 

8.  Computer  hardware  issues.  This  is  the  category  for  analysis  of  hardware  capabilities  and 
limitations  in  meeting  various  types  of  electromagnetics  computational  requirements.  Vector  and  parallel 
computational  techniques  and  implementation  are  of  particular  interest. 

Applications  of  interest  Include,  but  are  not  limited  to,  antennas  (and  their  electromagnetic 
environments),  networks,  static  fields,  radar  cross  section,  shielding,  radiation  hazards,  biological  effects, 
electromagnetic  pulse  (EMP),  electromagnetic  interference  (EMI),  electromagnetic  compatibility,  power 
transmission,  charge  transport,  dielectric  and  magnetic  materials,  microwave  components,  MMIC 
technology,  remote  sensing  and  geophysics,  communications  systems,  fiber  optics,  plasmas,  particle 
accelerators,  generators  and  motors,  electromagnetic  wave  propagation,  non-destructive  evaluation,  eddy 
currents,  and  Inverse  scattering. 

Techniques  of  Interest  include  frequency-domain  and  time-domain  techniques,  integial  equation 
and  differential  equation  techniques,  diffraction  theories,  physical  optics,  moment  methods,  finite 
differences  and  finite  element  techniques,  modal  expansions,  perturbation  methods,  and  hybrid  methods. 
This  list  is  not  exhaustive. 

A  unique  feature  of  the  Journal  is  the  publication  of  unsuccessful  efforts  in  applied  computational 
electromagnetics.  Publication  of  such  material  provides  a  means  to  discuss  problem  areas  in  electromag¬ 
netic  modeling.  Material  representing  an  unsuccessful  application  or  negative  results  in  computational 
electromagnetics  will  be  considered  for  publication  only  if  a  reasonable  expectation  of  success  (and  a 
reasonable  effort)  are  reflected.  Moreover,  such  material  must  represent  a  problem  area  of  potential 
Interest  to  the  ACES  membership. 

EDITORIAL  REVIEW 

In  order  to  ensure  an  appropriate  level  of  quality  control,  papers  are  refereed.  They  are  reviewed 
both  for  technical  correctness  and  for  adherence  to  the  listed  guidelines  regarding  Information  content. 
Authors  should  submit  the  initial  manuscript  In  draft  form  so  that  any  suggested  changes  can  be  made 
before  the  photo-ready  copy  is  prepared. 


The  ACES  Journal  is  flexible,  within  reason,  in  regard  to  style.  However,  certain  requirements  are 
in  effect; 

1.  The  paper  title  should  NOT  be  placed  on  a  separate  page.  TTae  title,  author(s),  abstract,  and 
(space  permitting)  beginning  of  the  paper  itself  should  all  be  on  the  first  page.  The  title,  author(s),  and 
author  affiliations  should  be  centered  (center-justified)  on  the  first  page. 

2.  An  abstract  is  REQUIRED.  The  abstract  should  state  the  computer  codes,  computational 
techniques,  and  applications  discussed  in  the  paper  (as  applicable)  and  should  otherwise  be  usable  by 
technical  abstracting  and  indexing  services. 

3.  Either  British  English  or  American  English  spellings  may  be  used,  provided  that  each  word 
is  spelled  consistently  throughout  the  paper. 

4.  Any  commonly-accepted  format  for  referencing  is  permitted,  provided  that  internal  consistency 
of  format  is  maintained.  As  a  guideline  for  authors  who  have  no  other  preference,  we  recommend  that 
references  be  given  by  author(s)  name  and  year  in  the  body  of  the  paper  (with  alphabetical  listing  of  all 
references  at  the  end  of  the  paper).  Titles  of  Journals,  monographs,  and  similar  publications  should  be  in 
boldface  or  italic  font  or  should  be  underlined.  Titles  of  papers  or  articles  should  be  in  quotation  marks. 
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